Automated Quantification of Hematopoietic Cell – Stromal Cell Interactions in Histological Images of Undecalcified Bone

Despite the progress in modern optical imaging methods, the analysis of histological data is still often hindered by the lack of proper quantification tools and methods, or by biased analyses that focus on a small area of interest. The synergistic approach presented here combines image analysis covering the entire bone marrow region, automated segmentation and object recognition of various hematopoietic and stromal cell types, co-localization analysis, and finally a validation tool of non-randomly occurring contacts provided by a new custom-designed simulation software.

Histology of whole bones including the bone marrow has been difficult to perform due to the various tissue densities present in one section – on the one hand the very hard, mineralized bone, on the other hand the extremely fragile marrow, which gets easily disrupted when cut together with the bone. As an alternative, the tape method for cutting bone sections presented here^14,15 has the advantage that it stabilizes the tissue, thereby leaving the endosteal areas intact (Figure 1). Besides being rich in osteoblasts, these areas have been suggested to play a role in stem cell maintenance²⁰.

The crucial role of stromal cells in developmental processes and cellular maintenance in the bone marrow is becoming increasingly evident. However, the analysis of these rare, fragile cells has proven to be difficult for two reasons: first, they are difficult to isolate from the bone marrow; second, the degree of heterogeneity among the stromal population is not known, and therefore one might miss an important population by using certain markers that have been described for stromal cells. The method to generate fluorescent bone marrow chimeras is useful to fluorescently mark and subsequently analyze the bone marrow stromal cell population as a whole. However, it should be noted that in these mice all mesenchymal cells of the bone marrow are visualized, including endothelial cells and osteoblasts. This needs to be taken into account if data from such chimeric mice is analyzed. Even though there is also a risk of hematopoietic cells of recipient origin surviving the irradiation, these cells typically only cover a very small area per field of view compared to the donor CD45⁺ cells, thus not affecting the histological analysis⁹. Clearly, this method should be combined in a subsequent step with a more specific detection of stromal (sub)populations. Currently, the analysis of these cells is mostly done by histology due to the reasons mentioned above, imposing a limit on the number of parameters that can be analyzed at the same time. The development of multi-parameter analyses in microscopy will help to overcome these limits.

During image analysis, the segmentation was performed using the original Otsu’s algorithm. Essentially, this clustering-based thresholding method automatically determines the intensity threshold for the optimal separation of the foreground (signal) and background (noise) pixels into two groups with the lowest statistical divergence in a given channel, resulting in a binary image¹⁷. Here, Otsu's method was applied to the natural logarithm version of the original image, that is:

Input image for Otsu = ln(Input image)

This works better for fluorescence images since Otsu's method may identify dim cells as background. Adding the logarithmic function allows for a better separation of the cells and background into two different classes by Otsu's algorithm. However, intensity thresholding alone is not enough to precisely detect single objects. It works efficiently when detecting well-separated objects in images, but it is not sufficient to separate single cells located within clusters. As some of the cell types that were of interest for us, such as eosinophils or B cells, are frequently clustered together in the tissue, their nuclei in the DAPI channel were segmented using an algorithm based on the k-means clustering of the intensity histogram with k = 10 ²¹. The cell segmentation algorithm will plot the brightest to the dimmest pixels (clustered in 10 bins) and constantly look for the formation of cell-resembling objects. These objects will then be defined as cells. When all pixels have been used, the resulting nuclei are dilated to generate a mask that is then applied to the other channels for cluster separation.

It should be noted that this analysis works well for hematopoietic cells that are compact, but that it is limited with respect to stromal cells, as it is not possible to determine the borders of these large, spread-out cells due to the interconnection of their dendritic extensions that form a reticular network. Thus, cell counts for stromal cells cannot be provided. Labeling single cells via fate-mapping reporters under the control of stroma-specific promoters will solve this issue, according to what has been published for neurons (using the “Brainbow” system²²) and for the labeling of follicular dendritic cells in the lymph nodes²³.

In addition to applying a noise reduction step to all channels by excluding objects below 30 px, size exclusion was applied in the process of object recognition for hematopoietic but not stromal cells. Small fragments of cells that were cut by sectioning may be excluded; however, the loss can be neglected in this case in favor of an unambiguous recognition of the cells. Extending the analysis and simulation from two to three dimensions will overcome this limitation. The development of methods to study cellular distribution in whole mounts, e.g., by multi photon microscopy or light sheet microscopy/SPIM, will certainly help in this respect. In order to study the complex cellular composition of multi-component niches, it will also become necessary to extend the number of parameters that are being analyzed in situ. Promising approaches for multiparameter analysis of histological information in combination with cytometric quantification have been published recently^24,25.

The contact analysis was performed by quantifying the overlap of two objects. Direct contact was defined as an overlap of at least one pixel. Importantly, this analysis is limited by the resolution of the microscopy images and therefore depends on the wavelength and mode of excitation, as well as on the numerical aperture of the objective lens used in the microscope, and the pinhole size. In the analysis shown here, a 0.8 NA, 20X objective with fluorochrome combinations excited with 488, 561, 594, and 633 nm using one photon excitation was applied. Hence, lateral resolution values between 0.372 and 0.483 µm were achieved. This allows assertions to be made about the juxtaposition of various cell subsets; however, it does not give any information about putative molecular mechanisms involved in cellular interactions. Additionally, other methods such as 2-photon intravital microscopy are needed to improve the characterization of these intercellular contacts⁹. Förster-resonant energy transfer (FRET)-based systems can help to confirm if these interactions are actually functional²⁶.

In order to evaluate the cellular microenvironment of tissue niches, it is not only crucial to analyze the direct cell contacts, but also to quantify the cell types in the neighborhood of a certain cell type of interest (“vicinity analysis”). Previously published studies measured distances between cells and tissue structures such as vessels based on the position of the nucleus²⁷. Since we were interested in determining the distance of various cell types, in part with varying ratios of nuclear to cytoplasmic volume, we preferred to measure the distance between surfaces. To this end, we calculated the vicinity radius by determining the Euclidian distance from a cell´s boundary after converting the values from µm to px. The Euclidean distance is the straight-line distance between two pixels and can be described by the formula²⁸:

Cells with pixels of their boundaries within an Euclidean distance of e.g., 10 µm from the pixels of a target cell’s boundaries were counted as neighboring this cell.

In its current form, the analysis only determines whether a cell is contacted or approached by another cell of the same or different type, thus providing only a binary YES/NO characterization. The actual number of cells that contact the cell of interest or are within a certain radius from it is not calculated. The next step will be to extend the current analysis in a way that the group size of contacting or neighboring cells can be detected and compared between different target cells of the same type. This will lead to important additional information about the composition of bone marrow niches, such as the survival niche for bone marrow plasma cells, for which a varying contribution of hematopoietic accessory cells has been discussed^29-33.

The object recognition algorithms were optimized together with the specialists from Wimasis using their commercially available custom solution service and are available upon request. Another option is to use either commercially available software for image analysis, such as Definiens, Volocity, or Imaris, or freeware such as CellProfiler.

The automated cell segmentation and image analysis tools were validated by comparing their results with manually analyzed data provided by 6 trained image raters with respect to two parameters, namely the object size and the number of cells in various cell populations. As shown in Figure 3B, the average object size for hematopoietic cells (PCs, eosinophils, and B cells) was similar between both methods, with a higher variance for manual size determination for plasma cells and eosinophils, probably due to their more irregular shape compared to B cells. The automated detection of cell numbers yielded slightly lower values than the manually detected ones. Notably, they were still within the range of the manual values in the case of eosinophils and PCs, while they were below these values in the case of B cells, probably due to a higher heterogeneity of expression of B220, which was used as a marker for B cells. B220 is known to be up-regulated during B cell differentiation in the bone marrow, and B cells with low as well as high intensity staining for B220 can be observed. B cells with low B220 signals fell below the applied threshold, which may have led to the exclusion of the earliest-stage B cells. A lower threshold for the automated detection may solve this issue. For stromal cells, which are harder to detect due to their less compact, more dendritic, and stretched-out morphology, the manual detection resulted in high variances between individual raters. Interestingly, the average total area of stroma was equal to the total area determined automatically, indicating that the analysis is well suited to compensate for individual bias by the raters. Taken together, these data show that the automated results generally represent the average values of manual analyses and are well suited to reduce interrater variability in the analysis.

In order to discriminate random and non-random positioning of cells within the structural limitations of the tissue, a simulation tool was developed. During the simulation, an artificial image with randomly positioned hematopoietic cells is created for every recorded bone marrow histological image. First, the stroma channel image is utilized in its original format. A mask is generated from the DAPI image by applying a fixed intensity-based threshold, thus converting the image into binary format; the binary mask then undergoes multiple steps of pixel-based dilation and erosion. The mask defines areas in the image field where simulated cells are allowed to be placed. The coordinates of the centers of the simulated cells are provided by an uniform random number generator, the boundaries of which where determined by the image size. Information about the cell number and average cell size for each population determined by automated image analysis of the recorded images are used to define the simulated hematopoietic cell populations. The simulated cell diameter is assumed to follow a one-dimensional Gaussian distribution from which the actual cell diameter is selected randomly: The diameter is then modified so that cell size cut-offs for the minimum allowable object size are introduced. In case the simulated cell would be positioned such that one or more pixels overlap with a no-go area or would be within a predefined distance from another already placed cell (4 µm from center to center as determined in recorded images), new random coordinates will be chosen whilst the diameter is left unaltered. This will be repeated until the number of simulated cells equals the number of cells in the recorded image.

The tool was based on the assumption that hematopoietic cells can be positioned freely within the bone marrow, while stromal structures act as a fixed matrix within the tissue. The tool was optimized to the situation in the bone marrow where parenchymal areas are crossed by a complex system of sinusoids. A similar modeling approach was recently published by Frenette and co-workers in order to analyze the positioning of hematopoietic stem cells in relation to stromal cells and vascular structures, which make up about 30% of the total femoral bone marrow volume²⁵. In order to correct for this effect, the simulation tool presented here is based on a mask of areas were cells are allowed to be placed. This was achieved by using a nuclear stain image, in which the objects representing the nuclei were expanded by a 5-step dilation (see Figure 6) after binarization. Areas with a low cellular density, i.e., low density of nuclei such as bone and sinusoids, will not contribute to the mask and therefore become no-go areas. In order to avoid an artificial reduction of these no-go areas by the dilation procedure, a 5-step erosion was applied subsequently, thus reopening larger exclusion areas, whilst leaving the high nuclear density (hence “allowed”) areas unaltered. This method can be easily adapted for other lymphoid organs. In those tissues were this method might not work (e.g., in areas where cells with a high cytoplasm/nucleus ratio are present), other methods such as cytoplasmic labeling can be used to create the mask.

Hematopoietic cells were simulated as circular objects whose size was determined by using a Gaussian random number generator, the mean of which agreed with the calculated average diameter of each cell type. The width of the distribution was based on measured cell size distributions in recorded images as determined by an image analysis software. In order to take into account that small cellular fragments were excluded in the automated image analysis, cell size cut-offs were introduced as determined from the measured cell size distributions for each cell type (see Figure 5).

Naturally, this works well for cell types that are relatively uniform in their size and shape but may lead to problems in the case of cells that are highly heterogeneous with respect to these parameters. It should be noted that this analysis is optimized for the major hematopoietic cell types in the bone marrow (granulocytes, hematopoietic progenitors and lymphocytes), which have in common that they possess a regular, nearly rounded shape. However, this method has not been tested for dendritic cells or macrophages, cell types with strongly irregular cell bodies. Analyzing such cells will present new challenges to our automated image analysis as well as to the simulation approach, including the need for improved cell cluster separation algorithms²⁴ as well as detailed shape-analysis and simulation of not only differently sized but also differently shaped cells. An alternative would be a simulation approach that works with the exact shapes as recorded. However, it must be noted that this method would significantly increase the deterministic behavior of the model system.

For the simulated experiments, 1,000 repetitions of the simulation for every recorded image were generated. Using this many repetitions will reduce the relative error contributed by the simulation process itself to about 3%. Namely, if the cell-cell co-localization values are calculated from simulations that are random, then the simulation process itself will be characterized by an N^-1/2 error function, where N is the number of simulations per measured image. This is the contribution by the simulation process itself to the accumulated error of the measured co-localization values. Thus, for N = 1,000 the contributed error is 3.16%. The simulations shown here ran for 20 to 60 min per image per 1,000 repetitions when only saving the simulated images as .rect files (a simple text format) rather than as actual images in .jpeg or .tiff format (also an option in the software).

Taken together, this approach allows for an unbiased, high-throughput analysis of histological images and enables the generation of statistically relevant data. This may be useful in comparing cellular localization under different conditions. Moreover, the temporal component can in the future be integrated into a modeling approach of bone marrow cell movements, as more data on the motility of various cell types in the bone marrow becomes available. This method can then be used to simulate perturbations of the system, such as the mobilization of cells from the bone marrow into the circulation, e.g., neutrophils in the case of an acute inflammation³⁴. To further elucidate the function of the bone marrow as a storage place for immune memory cells, one can also imagine extending it to model the competition for survival factors. Also, the possible crosstalk between distinct niches could be addressed, as well as the induction of niches. Together, this will help to understand the physiological changes (e.g., triggers for regenerative behavior in stem cells) as well as the pathological perturbations of the system as a whole, for example in the case of autoimmune disorders, neoplasia, or acute inflammation. As part of an iterative concept of experiment and modeling, new hypotheses can be generated based on the simulations, which in turn can again be tested experimentally, thereby helping to further understand the function and interaction of various cell types within the complexity of tissues.