Automated Cell Detection and Morphometry on Growth Plate Images of Mouse Bone

Microscopy imaging of mouse growth plates is extensively used in biology to understand the effect of specific molecules on various stages of normal bone development and on bone disease. Until now, such image analysis has been conducted by manual detection. In fact, when existing automated detection techniques were applied, morphological variations across the growth plate and heterogeneity of image background color, including the faint presence of cells (chondrocytes) located deeper in tissue away from the image’s plane of focus, and lack of cell-specific features, interfered with identification of cells. We propose the first method of automated detection and morphometry applicable to images of cells in the growth plate of long bone. Through ad hoc sequential application of the Retinex method, anisotropic diffusion and thresholding, our new cell detection algorithm (CDA) addresses these challenges on bright-field microscopy images of mouse growth plates. Five parameters, chosen by the user in respect of image characteristics, regulate our CDA. Our results demonstrate effectiveness of the proposed numerical method relative to manual methods. Our CDA confirms previously established results regarding chondrocytes’ number, area, orientation, height and shape of normal growth plates. Our CDA also confirms differences previously found between the genetic mutated mouse Smad1/5CKO and its control mouse on fluorescence images. The CDA aims to aid biomedical research by increasing efficiency and consistency of data collection regarding arrangement and characteristics of chondrocytes. Our results suggest that automated extraction of data from microscopy imaging of growth plates can assist in unlocking information on normal and pathological development, key to the underlying biological mechanisms of bone growth. Corresponding author. M.-G. Ascenzi et al.


Introduction
Microscopy imaging of mouse growth plates is extensively used to assess development and potential pathology.Such imaging confounds current automated methods for cell detection.Indeed, application to microscopy images of growth plates of each of the classic methods of image segmentation and processing (e.g.Canny segmentation [1], cartoon-texture decomposition [2], k-means clustering [3]) does not take into account that: 1) the color intensity of each stain used to visualize a specific biological component can vary throughout the growth plate; 2) the characteristics (e.g.size, shape, orientation) of cells (chondrocytes) vary greatly within healthy, normal growth plate; and 3) the shades of colors within cells are present outside cells.Moreover, the cells (chondrocytes) positioned beneath the plane of focus of the image appear faintly, yielding a highly non-homogeneous background to the cells on the plane of focus.Current cell detection algorithms apply to simpler images that show either a more homogeneous background or the presence of a cell-specific characteristic (e.g.intensity or nucleus) usable as a starting point for cell detection [4]- [9].
Analysis of chondrocytes' number, size, orientation, height and shape offer insights on the developmental effects of repressing various genes, the lack of which can lead to bone growth disorders (e.g.[10] [11]).Longitudinal growth of bones is the result of a process involving cell division, migration, and then ossification that occurs in growth plates located at both the proximal and distal ends of the long bone.In healthy growth plates, the chondrocytes sequentially differentiate into four sub-types, including resting, proliferating, prehypertrophic and hypertrophic chondrocytes, stratified from each end of the bone towards its center (for reviews, [12] [13]).Throughout the growth plate, only some nuclei are visible with appropriate stain.The resting zone contains relatively quiescent small, round, densely packed cells.Upon entry in the proliferating zone, cells elongate medial-laterally and undergo division by mitosis.At each division, the newly produced daughter cell remains closely situated with respect to the mother cell but may be far apart from other cells.The cells begin to form stacks.In pre-hypertrophic and hypertrophic zones, chondrocytes become arranged in columns.Such cells begin to enlarge and express Indian Hedgehog while the expression of Sox9 is reduced [14].In the hypertrophic zone, the terminal enlarged chondrocytes are larger than that in the rest of the growth plate, either round or elongated in the longitudinal direction, and packed closely to one another.The bottom of this region is marked by ossification [15], [16].Correlations exist among height of hypertrophic chondrocytes, growth plate length, and limb length [17].
In this paper, we propose a method of automated multi-step image processing.These steps prepare an image for automated measurement of the characteristics of the chondrocytes located on the plane of focus of the original growth plate specimen.Rather than manually determining cell profiles, automated cell detection, and subsequent automated morphometry would aid biological research by increasing efficiency of measurements.

Microscopy
The ideal image should be of a specimen whose preparation involves a stain that maximizes the contrast between cell cytoplasm and background, and in focus throughout.We visualized chondrocytes in the growth plate of mice ex vivo with a BX60 microscope (Olympus America) set to bright-field.The animal protocol was approved by UCLA's Animal Research Committee.Male mice on mixed background, predominantly 129/SvJ crossed with C57Bl6/J (WT for wild type), were sacrificed at perinatal 0, two weeks, four months.After isolation of hind limbs, we stained coronal sections of femora and tibiae to visualize cartilage (with Alcian blue or Safranin O), cell cytoplasm (Fast Green) and nucleus (Hematoxylin) [18] [19].We took TIFF images of sections at 60× and 100× magnification.

Image Preparation
We imported the TIFF images of growth plates into software XaraX version 1.1 (XaraX Group).The Shape Edi-tor, Combine Shapes and Intersect Shapes allow tracing of the region of interest (ROI), cropping along the traced line to eliminate the image of the tissue surrounding the growth plate and possibly select a region of the growth plate, and replacement with a white background (Figure 1).This step is necessary to avoid interference of colors that do not belong to the ROI during algorithm application.We then exported the cropped image as TIFF image with same resolution as the original image.

Cell Detection Algorithm (CDA)
Our CDA is written with and runs on MATLAB version R2012a (Math Works), and is available for download as file CDAlg.m together with three of the full resolution images on which the CDA was run to prepare some of the figures in this paper, from https://sites.google.com/site/bthsca/automated-cell-detection-1.The input for the CDA consists of a TIFF image (Figure 1).The output is a black and white (i.e.binary) image of the size of the input image, with white pixels representing cells on a black background.The CDA consists of 5 steps: Step 1: We assign similar intensity values to cells in focus while we render the background, defined as cells out of focus and extracellular matrix, more homogeneous.We use the color blending method of Retinex, originally developed to simulate the visual cortex blending of colors [20] [21], and then widely applied to imaging [22] [23].We use the formalization of Retinex in terms of a Poisson equation [24].Given an initial image, f, we seek a reconstructed image u with is the discrete Laplacian at pixel, (i, j) and where is the hard thresholding function and τ the thresholding parameter.Thus, if a small gradient (of magnitude smaller than τ) is present at a pixel, it is replaced with zero.The Poisson equation then constructs the image whose gradient most closely matches the vector field that models the difference between each pixel and its neighboring pixels.We use the following semi-implicit scheme that converges rapidly to the steady state for local variations: ( ) 4 We apply this process iteratively, NI (number of iterations) times, each time to the previous resultant image.
Step 2: We enhance cell profiles (edges) with anisotropic diffusion.The anisotropic diffusion function needs to discriminate intensity changes by location: near vs. far from edges, and along vs. across edges [25].Let f be an image, u the filtered image and D the diffusivity tensor.We used the numeric approximation for single time-steps with g ( ) and 1 set as eigenvalues of D.
Step 3: We extract cell profile information with gradient thresholding.We use the gradient matrix generated in Step 2, containing values scaled between 0 and 255.The hard threshold parameter, called S for "separation", is chosen to identify single cells within packed clusters, while minimizing pixel loss from cell profile.We use the same Equation ( 4) with S replacing τ.We obtain a binary representation of the image with white cell profiles on black background suitable for subsequent morphometric analysis.
Step 4: We use the convex hull operation [27] to identify the smallest convex set of points that fully encloses the profile of each cell.This operation renders the profile connected and filled because the thresholding needed to obtain the binary representation at Step 3 may have caused interruptions of the cell profile.This operation is appropriate because 1) cells show an elliptical and therefore convex shape [28]; and 2) the number of missing pixels from the profile, minimized in Step 3, is small.However, when two cells in focus are close to each other, the two separate cells may lump together and their profiles may form a larger concave shape.
Step 5: We eliminate profiles that do not reflect actual cell size or shape with hard tresholding.Overly large or overly small or misshaped profiles are often created during convexification because of proximity of incomplete profiles or proximity of an incomplete profile with the boundary of the growth plate region.Alo (area lower) and Aup (area upper) bound in pix 2 define a range out of which the profiles are excluded.Because these boundaries depend on image size and magnification, they need to be set by the user.Further, we use set up shape thresholding in terms of isoperimetric ratio (IR) where A denotes area and P perimeter of the shape.The threshold was set at 0.5 because an NI greater than 0.5 demonstrate geometriesS different from ellipses.
The proposed CDA, implemented with MATLAB R2012a with 64-bit Windows-7, an Intel @ Core TM i3-370M processor (2.40 GHz) and 4 GB of RAM, takes between 5 3 10 − × and 5 7 10 − × sec/pix to create the output image with memory usage of 0.000149 mb/pix (Appendix 1).The computational cost O(MN) of running the CDA on an M × N image is at most equal to HMN for some finite constant H, as M and N approach infinity (Equation ( 9) of Appendix 2).

Manual Cell Detection (MCD)
MCD was conducted on 390 cells on a total of three images with XaraX software, for comparison with the CDA.Because chondrocytic two-dimensional profiles throughout the growth plate are best modeled as ellipses [28], we used Ellipse and Rotation tools to fit a white ellipse to each cell profile that appeared in focus.Because we traced only cells in focus, there is no issue of tracing cells that appear to overlap.By drawing a black background as large as the original image, we obtained a black and white image that was exported at the original image resolution.

Automated Morphometry of Images
We conducted morphometry on images from each of CDA and MCD with MetaMorph software (Molecular Devices).The Color Threshold was set from 127 to 255 for each of Red, Green and Blue.We calibrated Meta-Morph in conjunction with the magnification used to generate the image to measure lengths in real microns.The output parameters were cell number, area, orientation (measured with respect to the long bone axis, counterclockwise from −90˚ to 90˚), height, and shape in terms of isoperimetric ratio IR (see Equation ( 8)).

Statistical Analysis
The robustness of the morphometric MCD was measured in terms of inter-and intra-observer errors.The magnitude of these errors was assessed for cell number, area, orientation, height and IR.To assess the magnitude of such errors, the cell profile of six randomly chosen cells in focus was manually drawn with XaraX and measured with MetaMorph forty-two times by each of two observers, twice, two months apart.Further, to detect significant differences between CDA and MCD, and between adjacent zones of a given image, we used the t-test for each of the measured parameters, and the Dice index to measure similarity of cell shapes [29] [30].Because the sample size is greater than 40, the possible skewness of the data distribution does not affect the robustness of the t-test.We consider a p-value smaller than 0.05 as statistically significant.The collected data are presented as mean ± standard deviation.

Results
We have developed an automated method for analysis of growth plate images.Figure 2 shows intermediate images and output image of the CDA.

CDA
The CDA depends on the five parameters τ, NI, S, Alo, and Aup whose values depend on the image characteristics: τ (Step 1, Equation ( 4)) is the threshold value at which small visual differences between a given pixel and its neighborhood are eliminated.At small values, we obtain numerous white dots corresponding to either miscellaneous background details or cell fragments.At high values, all cells are homogenized into the background and therefore not detected.NI, the number of iterations of Step 1, has an upper bound, beyond which the cell profiles become too blurry for correct detection.When NI is appropriately smaller than such upper bound, it sufficiently homogenizes the background noise without blurring the cell profile.S is used to identify the cell profile (Step 3).S has an upper bound below which it identifies cell profiles with small pixel gradients.The use of small pixel gradients has the benefit of identifying the complete cell, at the expense of multiple cell aggregation.Higher values of S tend to identify only the high contrast between the nucleus and the stained background leading to reduced cell aggregation, but also to partial cell detection.For stains that do not provide good contrast between cell on plane of focus and background, a low value for S is necessary to detect the small pixel gradients defining the cell profiles.For stains that provide high contrast between the majority of the cells in focus and background, a higher value of S is suitable.Alo and Aup (Step 5) control the area of the detected shape: Alo eliminates leftover minute fragments and Aup eliminates oversize shapes that cannot represent cells.
With appropriate values of the parameters, the CDA analysis confirmed known facts regarding the mouse growth plate in tibia and femur [16] [31].For the WT, the mean cell area is significantly smaller in the resting zone than in the proliferative zone.Table 1 shows data from a femur as mean ± standard deviation.Significant differences are shown between adjacent zones of the growth plate with ^ for p < 0.01 The values of Alo was optimized for the resting zone that contains the smallest cells.The large standard deviations of cell area are due to proliferation that significantly increases cell size and to hypertrophy that further increases cell size.The significant differences in mean cell area among zones results from distributions of resting, proliferative and hypertrofic cell areas with tails that become significantly longer as cell size increases from resting to hypertrophic zone.This is also the reason why cell area from the three zones ranges shows only a small partial overlap, and give rise to significant differences in cell areas between adjacent zones.
The CDA captures more cell area (per cell) in the proliferative zone because it picks up the parts of the cell which are stained: the CDA is far more sensitive to small gradients in color change, hence it is better able to pick up the color transition between the background and cytoplasm even if the stain bleeds into the cell.We  found that the mean cell area is significantly smaller in the proliferative zone than in the hypertrophic zone; and that the mean orientation differs significantly between resting and proliferative, but not between proliferative and hypertrophic zones.Further, the mean isoperimetric ratio IR is significantly different between resting and proliferative zones, and between proliferative and hypertrophic zones, indicating that the cells become less circular as they progress from the resting to proliferative, and from proliferative to hypertrophic zone.Furthermore, our CDA confirms the rounder shape, and the orientation off the horizontal (medial-lateral) direction, of chondrocytes in the Smad1/5 CKO mutant in comparison with the control WT mouse (Figure 3), previously determined by manual detection [10].Note that for the CDA analysis of this WT growth plate, the Retinex parameter NI was set to zero because these fluorescent microscopy images showed no faint cells behind the plane of focus.In fact, Retinex color blending on inherently high contrast fluorescent microscopy images would reduce the contrast between cell profile and background, making cell detection sub-optimal. .Orientations are marked in (e, f) to emphasize differences.

MCD
The 390 cells were chosen on the plane of focus of a total of three images by human eyesight judgment.The chosen cells showed presence of the 92 to 100 range of Red, 62 to 100 of Green and 55 to 100 for Blue.The out-of focus cells appearedoverall darker with 86 to 100 range of Red, 49 to 97 of Green, and 45 to 100 of Blue.The inter-and intra-observer errors for a unique measurement were found to be smaller than 1.5%.We show the data of the calculation of inter-and intra-observer errors concerning the cell area from Figure 1 (Table 2).The minimum (min), maximum (max), mean, and standard deviation (stdev) of distributions of 42 measurements of 6 cell areas provide data for calculation of comparison relative to observer (Obs) 1 and 2 during first (I) and second (II) setting, in terms of % errors (in italics) for measurements considered either paired (p) or unpaired (u) All % errors are smaller than 1.5%.The error of a unique measurement is at most equal to either the difference of the two largest measurements divided by the smallest measurement, if the measurements of the two sets are considered paired (comparison of two measurements of the same cell), or to the largest measurements' difference divided by the smallest length, if the measurements of the two sets are considered unpaired.On the basis of the small magnitude of the above-analyzed morphometric errors, performance of one measurement by a single observer was deemed appropriate.The computation of errors on six cells per image remained smaller than 1.5%.

Comparison between CDA and MCD
The best detection of the CDA in comparison to MCD occurs with a 5% difference on 390 cells detected on three images, a Dice similarity index of 0.88 in cell number, and no significant differences for cell position, area, orientation, height, shape factor in terms of isoperimetric ratio.In particular, this means that the convexification step in the CDA did not produce a shape with significant distortion.The aggregation of adjacent cells was present in all zones of the growth plate when single values were assigned to each of τ and S for the whole growth plate.To decrease such aggregation, the parameters τ and S need to vary across the growth plate, specified for each of resting, proliferative and hypertrophic zones, due to variation of relative distance among cells.The hypertrophic zone presents the highest challenge to the CDA because the cells are tightly packed and their profiles frequently touch (Figure 4 [17]), we found small variations in hypertrophic cell height and growth plate length to mean hypertrophic cell height ratio (Figures 4(c)-(e); Table 3) for different choices of CDA parameters.In  , for which the best cell detection occurs (shown in italics).For all runs, NI = 85, Aup = 20,000.Figure 5 shows another example of CDA running with different parameters' values that focuses on the middle of a primary ossification center (bottom) and of a secondary ossification center (top).

Discussion
We have developed a specialized algorithm, implemented in MATLAB, which performs the automated detection of cells in growth plate images.We have used the algorithm to verify known biological characteristics of the mouse growth plate and provide a tool for further research regarding, e.g., implications of height of hypertrophic chondrocytes on limb length.
During the process of developing CDA, we considered current methods for image detection.Canny segmentation, cartoon-texture decomposition (Figures 6(a)-(c); [1] [2] [32]- [34]) created false positive cell detection from intensity variations of stain.Intensity variations produced also inconsistent clustering of background pixels and cell pixels with k-means clustering (Figure 6(d); [2] [3]), which clusters data points on similarity basis.Current cell detection algorithms, such as Cell Profiler [4] and HT-HCS [5], do not apply to our images because they are based on existence of "primary objects" (such as presence of cell-specific intensity), one per cell.We quantitatively compared our CDA to image enhancing processor Image J [6] that includes various filtering, thresholding, segmentation compared our CDA to image enhancing processors Image J [6] that includes various filtering, thresholding, segmentation and clustering methods [7], with a recent implementation to detect single cells on a slide [8].We compared each of Image J and CDA ( ) to MCD of 390 cells from three images.We obtained the best readable image with Image J, by first rendering the background more homogeneous with a 50% contrast enhancement.Then to obtain a black and white image, we used color threshold "HSB color space" to threshold black and white.Afterwards, we adjusted hue, saturation and brightness to maximize individual cell detection.We found that Image J picks up background cells out of focus and often separates visible nuclei from cells.Therefore, Image J generates at least a 26% difference in cell count in comparison with the MCD results.This percentage is larger than the 5% difference generated by the CDA with respect to the MCD.Image J showed no significant difference with MCD for cell area, orientation and height while the shape factor (isoparametric ratio) differed significantly between Image J and MCD.This in comparison to no significant differences for cell area, orientation, height and shape factor on all three images obtained with the CDA.Further, any CDA output image can be processed with existent software to assess more complex parameters for object and pattern assessment [35]- [38].
Recently, Buggenthin et al. [9] developed an algorithm for detection of hematopoietic stem cell in culture on  a gray-scale image obtained under bright field microscopy.This algorithm uses a machine-learning based background correction method, which plays a role analogous to that of Retinex and anistropic diffusion in our CDA.After correction, both CDA and Buggenthin algorithm use a hard-thresholding method to identify cells.The Buggenthin algorithm includes a marker based water-shedding method to separate the clumped cells, while our CDA lacks such step.Therefore, when our CDA is applied to the image of hematopoietic stem cells, aggregation of cells occurs, and 5% of cells, small ones, were not detected.However, when we applied the Buggenthin algorithm to our colored growth plate images, it does not perform as well as the CDA.In fact, it produces a single-colored image of pure background because our colored growth plate images have a considerably high degree of background heterogeneity, and the color of many cells is not significantly different from their neighborhood background.
Complexity and local variations of the images are such that the automation of the CDA does not extend to the choice of the values of the parameters τ, NI, S, Alo and Aup that need to be chosen by the user.In fact, the optimal values of these parameters can only be determined by user's experiment and experience on the specific image.We note that the cells on focus need to be resolved enough on the image to be clearly visible by human eye against the background, for the CDA to be able to detect them.Presence or absence of stain of the nucleus does not interfere with appropriate cell detection with the CDA.In fact, while the black stained nucleus causes detection of only the cell around the nucleus, the convex hull operation fills up the space occupied by the nucleus, therefore providing results similar to the results obtained for images with unstained nucleus.Also, if the  CDA leaves a few cells lumped together, a quick analysis of the morphometry data outliers can assess whether such lumping affects the mean height of the distribution.Such outliers can be excluded if needed.

Conclusion
In conclusion, we have presented here an innovative algorithm that performs accurate cell detection of growth plate images.In particular, our CDA can provide rapid measurement of height of hypertrophic chondrocytes, for comparison to height of hypertrophic zone, of growth plate size and limb length, for a first-screening before investigations on phases of volume increase.This method is applicable to specific genetic mice models to investigate the biological mechanism of limb lengthening. of size M × N. The computational cost is MN.
10) Divide the y from 9) by 2. The computational cost is MN.
Step 3 has a computational cost of MN because the gradient threshold compares each pixel's color code gradient to the gradient threshold.Step 4 has a computational cost of P(MN) because the number of the convex hull operations is at most MN.The black pixels that are inside the boundaries of each polygon generated through convex hull operation are then transformed into white pixels through MATLAB built-in function "poly 2 mask" and the computational cost is P(MN).
Step 5 has also a computational cost of P(MN), since the number of polygons examined for size and shape thresholds equals the number of convex hull operations.

Figure 1 .
Figure 1.Image of two-week old WT growth plate viewed with bright-field microscopy after excluding surrounding tissue (inset).Cytoplasm is stained gray-green with Fast Green; nuclei are black with Hematoxylin; cartilage, mucin and mast cell granules are pink-red with Safranin O; and bone is blue with Fast Green.

Figure 2 .
Figure 2. Input, intermediate steps and output of the CDA obtained with 7, 85, 130, 160 NI S Alo τ = = = = and 20,000 Aup = .(a) Original image, detail of Figure 1; (b) Image resulting from reduction of visibility of cells out of focus, and of inconsistencies in background shading and luminosity (Step 1); (c) Cartoon-like component image resulting from application of low diffusion across edges and high diffusion at locations of noise and small features (Step 2); (d) Output image resulting from gradient thresholding to address cell separation (Step 3), convexification (Step 4), and size and shape thresholding (Step 5).Cell clustering can be reduced or eliminated by choosing zone-specific values for and S (see Figure 4).
(a), Figure 4(b)).The analogy of best images for different values of τ and S. (Figures4(c)-(e)) is indicative of the synergistic effect of τ and S and that the parameters' values are not 1-to-−1.Within the context of the relation between hypertrophic cell height and length of growth plate, bone and limb ([4]

Figure 4 .
Figure 4. Specification of parameters τ and S for hypertrophic zone of tibia's growth plate of neonatal (P0) WT mouse.(a) This is the original image of growth plate where the cartilage is stained with Alcian blue and the nuclei are stained black with Hematoxylin.Because a specific stain is not used for the cytoplasm, the blue of the cytoplasm is due to the Alcian blue stain.The resting zone is delimited by red and the proliferative zone by green.The hypertrophic zone is delimited by blue and enlarged in (b).We show CDA images obtained with (c) 7, 130 S τ = = ; (d) 5, 120 S τ = = ; (e) 5, 145 S τ = = .(c, d, e) yield detection of single cells within a 5% error of MCD.For all runs, 85, 400 NI Alo = = and 20,000 Aup = .

Figure 5 .
Figure 5. Specification of parameters τ and S for tibia growth plate of 4-month-old mouse.(a) This specimen image focuses on the middle of a primary ossification center (at bottom of image) and of a secondary ossification center (at top of image).The cytoplasm is stained gray-green with Fast Green, the nuclei are stained black with Hematoxylin, and cartilage is stained pink with Safranin O.The cartilage is pink, not red, because of over-staining with Hematoxylin.We focused on region (b).The CDA results are shown (c) with 80 S = ; (d) 100 S = ; and (e) 120 S = .For all runs, 7 τ = , 85 NI = , 200 Alo = and 20,000 Aup = .Image (d) provides the best result for detection of individual cells of correct size.

Figure 6 .
Figure 6.Results from standard segmentation and processing algorithms.(a, b) Canny segmentation (with edge function with Canny option in MATLAB) applied directly to an input growth plate image results in false positive cell-detection and incomplete edges; (c) The cartoon component of the cartoon-texture decomposition (within the Bregman Cookbook in MATLAB), results in false positive cell-detection from cells out of focus and changes of intensity in the background; (d) Results of k-means clustering (built-in function in MATLAB) show unwanted clustering of background pixels.

Table 2 .
Assessment of inter-and intra-observer errors.

Table 3 ,
the data are presented as mean ± standard deviation.
* indicates significant difference in

Table 3 .
Cell count and height of hypertrophic chondrocytes from Figures 4(c)-(e).