The availability and usage of advanced fluorescence imaging techniques such as confocal and multi-photon microscopy have dramatically increased and enabled researchers to investigate biological processes with a high degree of spatial and temporal resolution. This includes, but is not limited to, studies ranging from static detection of the subcellular localization of proteins to dynamic tracking of fluorescent probes in single cells to intravital imaging of cell behavior in complex tissues and organs. These developments have been especially fruitful for fields like immunology, cancer research, and neuroscience where the system behavior is largely governed by dynamic cell interactions, for instance, by the interactions between lymphocytes and antigen-presenting cells in lymph nodes
1–3, lymphocyte recruitment to and interaction with tumors
4, and synapse-glia dynamics in the brain
5. Quantitative analysis of cell interaction behavior can considerably increase the information we can gain about the molecular mechanisms governing cellular communication processes or may be used to assess and quantify the efficacy of drugs. However, the development of computational high-throughput methods for automated and standardized quantitative analysis of the resulting 3-D image data has lagged behind the experimental advances
6.
Here, we describe a protocol for the quantitative investigation of cell-cell, cell-tissue or cell-pathogen interactions that uses a fully automated, high-throughput image analysis method. The approach involves four steps that are performed automatically without user intervention: First, the actual (true-positive) signal is separated from background (false-positive image elements) in each fluorescent image channel. Second, individual image objects are detected in the output data of the first step, which permits acquiring object numbers and size statistics and allows for removal of image artifacts based on prior information about, for instance, minimum cell size. The third step merges the different fluorescent color channels to obtain an unambiguous segmentation. This step is especially important for interaction analyses because interfacing objects of different type (i. e., different color) usually show spatial image overlap that an accurate interface analysis must account for. Finally, in the last step, the interface areas are computed.
Advantages and Disadvantages of our method
Previous approaches to interaction image analysis relied on semi-quantitative estimations using manual measurements of such features. For example, in conventional bone-osteoclast interaction analysis, data sets are sent to commercial labs and processed by personnel manually delineating cellular boundaries and interfaces
7–11. In contrast, the method presented here enables investigators to perform rapid and standardized analyses that do not require operator interpretation in the vast majority of cases and are especially suited to the quantification of differences between experimental groups in large 3-D data sets.
To assess the limitations and the usability of our protocol for different data types and image features we tested different types of typical confocal/two-photon image data (see details below) acquired by different experimenters using different microscopy platforms. Moreover, we performed sensitivity analyses by assessing the impact of the variation of the threshold selection parameter and the introduction of artifical image noise on the interface analysis. These tests showed that our approach works equally well for all data types we tested and that parameter/quality variations did not significantly impact the interface analysis results. However, we tested typical image data and noise. A limitation of the protocol (and segmentation approaches in general) is that excessive noise and/or high background signal (getting close to the intensity of the object signal) may lead to a failure of the ability of the method to separate noise/background from actual signal which subsequently generates incorrect interface quantification results. As a rule of thumb the protocol can be applied to image data whose intensity histogram meet the shape characteristics we describe. Vice versa, the quality of the analysis can be expected to deteriorate if e. g. the peak in the lower intensity region of the histogram representing the background `collapses' into the higher intensity `signal' region. However, we would like to point out that the main application of our protocol is the comparison of different experimental groups (see following section for examples) and that for this purpose (small) segmentation errors are acceptable as long as they occur consistently which is guaranteed by our approach.
Example applications
The software we use in this protocol was initially developed for quantifying osteoclast-bone interactions in order to analyze how a S1P
1 (sphingosine-1-phosphate)-receptor knock-out mutant and the immunomodulatory drug FTY720 (Fingolimod, an S1P
1 agonist) that affect a lipid phosphate chemosensing pathway alter osteoclast activity in situ and thereby affect bone homoeostasis
12. We present this as our main example application and then demonstrate the versatility of the approach by describing its application to two other types of data. In the second example we show how our approach can be used to analyze image data to quantify the extent to which dendritic cells decorate the fibroblastic reticular cell (FRC) networks on which T lymphocytes migrate within lymph nodes. The last example shows the applicability of the protocol for the investigation of interaction phenomena on a smaller scale: We analyze image data from recently published experiments
13 about the influence of a mouse T cell signaling mutation in the gene encoding SAP (signaling lymphocyte activation molecule (SLAM) associated protein, the cause of X-linked lymphoproliferative syndrome (XLP)) on cell membrane contacts of T and B lymphocytes (this SAP mutation is known to affect the development of humoral immunity by influencing the stability of interactions between these T and B cells).
All these examples have in common that their analysis requires a quantitative assessment of either the size and/or the spatial localization of interfaces between cells or spatial regions with certain properties relevant for the biological question at hand but differ in terms of size and shape of the cells and tissue structures.
The applications demonstrate the advantages of using an automated and standardized analysis method to obtain unbiased comparisons of data sets obtained under different experimental conditions. As illustrated by the images in and , background fluorescence and the overlap between different color channels can make manual segmentation a challenging and subjective task. Under such conditions, variations in background signal/noise levels between data sets or between color channels in a single data set can easily mask subtle but significant biological differences or artificially create the impression of such differences. The potentially clinically relevant results obtained for the quantification of osteoclast-bone interaction robustly demonstrated the ability of the immunomodulatory drug FTY720 to reduce the adhesion of osteoclasts to bone surfaces based on multiple data sets displaying typical variations in image quality. This effect was not evident from visual inspection of the microscopy images and could be rigorously quantified only because of the adaptive standardization our protocol offers.
Background to the image analysis program
The image analysis program performs the following steps automatically:
First module: Signal-Background separation by automated adaptive threshold segmentation The critical step in image segmentation is the separation of signal from background, which includes discrimination of different signals (colors) from one another when considering multicolor fluorescence datasets. The first module of the image processing software that is used here performs an automated adaptive intensity threshold segmentation based on the characteristics of fluorescence confocal or two photon image data. Such intensity histograms show a peak in the lower intensity region that represents the background signal, while the actual image signal is – relative to the form of the background part – approximately uniformly distributed over the (higher) intensity spectrum. We found this to be characteristic of all two photon and confocal fluorescence microscopy images we analyzed and distinct from most conventional (non-fluorescent) image data (as in
14) (). Because the exact shape and location of the peak and constant region vary depending on the image acquisition features, a stable and reliable segmentation method capable of yielding comparable results has to adapt to these variations that may occur between different as well as within individual data sets. The segmentation algorithm used here accomplishes this by defining the transition point between background and actual signal (i.e., the threshold) as the intensity value at which the curvature of the histogram graph is maximal (see section “Automatic adaptive thresholding: Maximum curvature estimation” for details). In contrast to fully (3-D)-global thresholding methods with single constant or user-selected cut-offs for the whole 3-D image data set, our approach computes the threshold for each image z-slice separately and is therefore capable of compensating for intensity inhomogeneities (or “attenuation”) along the z axis (for instance, in deep-tissue 2-photon imaging
in vivo).
Automatic adaptive thresholding: Maximum curvature estimation In the case of photographic images histograms often have a bimodal shape, with one peak representing background and the other signal portion of the image (
Supplementary Figure 1). In such cases it is obvious that the threshold selection criterion is to find the minimum between the two peaks. There is no such criterion that is similarly obvious and also physically plausible in case of two-photon/confocal microscopy image data and the transition between background and signal cannot be defined as a precise location similar to that in case of bimodal histograms. It can only be defined with some remaining degree of fuzziness in the region of the transition between the strongly decaying part of the peak in the lower intensity spectrum and the region with the relatively flat slope. We chose the maximum curvature and its approximation by using a derivative criterion for normalized histograms as threshold selection criterion because it represents a mathematically definable, reasonable estimate of the transition point at which the strongly negative slope of the background peak changes into the relatively flat slope associated with the signal region and moreoever, because it is computationally efficient (alternative threshold selection criteria are discussed in
14).
To compute the curvature of the (discrete) histogram in the relevant region a Gaussian model may be used to fit the histogram in the relevant region, i. e., between global maximum (or local maximum with highest intensity value in the background region) and intensity 254. The maximum intensity histogram bin 255 is omitted in all calculations because it is significantly higher than the other “true positive” intensity values. This is due to the aggregation of all intensities ≥ 255 in this bin. We also omit the 0 bin, because in the case of image data of a low overall intensity (and thus many 0-intensity voxels) its inclusion may interfere with the generation of a normalized histogram that exhibits the shape features required for threshold computation.
With the angle α,
f (
x) the histogram fit function and the line element
ds given by
the curvature of the graph of
f (
x) is given by
(). In our implementation we approximated this process by using the (discretized) derivative of the histogram
h(i): To determine the intensity threshold used for the segmentation, the derivative Δ
h/Δ
i is computed stepwise starting at
i=254 with decreasing
i until the transition from the constant image to the non-constant background region is encountered. Prior to the threshold determination, the histogram is smoothed using a Gaussian filter (0 bin is the maximum bin of intensity histogram) scaled by the histogram integral and multiplied by a factor of 255 in order to normalize for different cell/bone tissue coverage of the image area. Because of the common histogram shape features present in all data sets, this scaling/normalization procedure resulted in the value of the discrete derivative
to be approximately 1 at the maximum curvature point (). Our sensitivity analysis of the precise location of the selected threshold (see section “Robustness assessment: Impact of variations of the slope constant on the interface area results” for details) shows that this approximation is valid.
Supplementary Figure 2 shows an exemplary threshold segmentation result overlaid with the original image.
Second module: Object detection using connected component analysis Following the automated thresholding procedure, our approach uses a 3-D connected component labeling algorithm (see
15 and references within) to identify individual image objects. The connected component labeling permits computation of the numbers of objects (such as cells) and the corresponding object features such as volume, surface area, and location. It is also required for obtaining detailed information on cell-cell interaction behavior (e.g., the number of cells of type A attached to cells of type B) and for comparing features of potential interest of interacting vs. non-interacting cells. In addition to the utility of the connected component analysis for the quantification of object features, our software also uses this analysis to improve image quality by applying a threshold filter for cell volumes, thus allowing removal from the processed dataset of small pieces of cell debris or bright specks resulting from artifactual dye labeling.
Third and fourth modules: Channel merging and interface area computation using normalized intensity comparison and voxel-neighborhood analysis The third step of the analysis pipeline finalizes the segmentation by adaptively merging the previously separate segmentations of individual color channels. This is a pivotal step because most current image data involve objects (cells, pathogens, tissue components) that are labeled with different colors and whose fluorescent signals overlap when those objects are in close proximity, resulting in overlap in the segmentations of the individual color channels. Disregarding this issue may lead to distorted interaction results and even to the accidental removal of small objects, especially if one channel is significantly dimmer than the other. Our channel merging approach accounts for spatial overlap between neighboring or interacting image objects by combining the single channel thresholding results with an adaptive intensity comparison of the original image data based on normalization of each channel with its intensity profile. After the initial threshold segmentation (performed for each channel separately) each voxel that has been assigned both class labels (`overlap voxel') is evaluated by comparing the original image intensities in each channel normalized with the corresponding maximum intensity values. The voxel is then assigned to the final segmentation class that represents the channel with the higher relative intensity at that particular voxel. This approach therefore allows for an accurate estimate of the actual interfaces (see and
Supp. Fig. 3).
Finally, the interface areas are computed by identifying those voxels in channel one that are direct neighbors of voxels in channel two. We use the sum of the surface voxels as an estimate of the total contact area. This is sufficient because the analysis result is provided in terms of ratios of interface areas to total cell or tissue surface areas or in terms of comparisons of such areas between different experimental groups.
Robustness assessment: Impact of variations of the slope constant and noise on the interface area results The robustness of an image segmentation and analysis method towards variations of pre-set parameters and image quality is an important measure of its utility. Therefore, we performed a sensitivity analysis to evaluate the impact of variations in the slope value used to identify the boundary between background and true image data. This analysis showed that the average ratio of the interface area / total surface area in two groups (bone-osteoclast gene knock-out vs. wild-type data) with 20 data sets each (each data set consisting of ~10 x–y images of 512×512 pixels) changes only about 10% even when varying the slope value by a factor of 2 in both directions for both channels. More importantly, the differences of the average area ratios between the control and knock-out groups – which is the relevant result in this kind of analysis – show a negligible dependence on the variation of the characteristic slope (change of difference 1.1%; see ). These results show that the particular choice of the slope value is less important than using the same constant throughout the data analysis, which subsequently generates the proper adaptive threshold value. This ensures that variations of image quality within and between data sets are normalized, thus permitting the consistent comparison of different data sets.
In addition to testing our approach with different types of typical confocal and two-photon microscopy data we assessed the robustness of the analysis under the influence of a variety of noise levels. We generated different levels of artifical noise for an exemplary data set (image noise follows a Poisson distribution
16). The comparison of the analysis results yielded interface area differences of ≤5% between different typical Poisson noise levels (see
Supplementary Figure 4 for details). However, due to the combination of noise and background signal and the diversity of possible features of fluorescence microscopy data we cannot provide straightforward rules to determine a priori which data our protocol will work with. Low data quality may lead to significant misclassification (see
Supplementary Figure 5) and we therefore advise the user to perform tests and examine test results to estimate the suitability of the automated analysis using this protocol and/or to optimize the user-definable parameters for the particular data. However, for the main application of this protocol, the comparative analysis of interaction behavior, it is above all important that data processing is standardized and that the data sets that are to be compared are of roughly similar quality. In that case segmentation errors within reasonable limits are not problematic for comparative analyses using our protocol because they occur consistently.
The publisher's final edited version of this article is available at 
