Scientific Reports

Sci Rep. 7: 43693

Multiscale image analysis reveals structural heterogeneity of the cell microenvironment in homotypic spheroids

Alexander Schmitz¹*, Sabine C. Fischera¹*, Christian Mattheyer¹, Francesco Pampaloni¹, Ernst H. K. Stelzer¹

1. Physical Biology/Physikalische Biologie (IZN, FB 15), Buchmann Institute for Molecular Life Sciences (BMLS), Cluster of Excellence Frankfurt – Macromolecular Complexes (CEF – MC), Goethe Universität – Frankfurt am Main (Campus Riedberg), Max-von-Laue-Straße 15 – D-60348 Frankfurt am Main, Germany

a. sabine.fischer@physikalischebiologie.de

*. These authors contributed equally to this work.

Copyright © 2017, The Author(s)

DOI: 10.1038/srep43693

Published online: 03 March 2017

Abstract

Three-dimensional multicellular aggregates such as spheroids provide reliable in vitro substitutes for tissues. Quantitative characterization of spheroids at the cellular level is fundamental. We present the first pipeline that provides three-dimensional, high-quality images of intact spheroids at cellular resolution and a comprehensive image analysis that completes traditional image segmentation by algorithms from other fields. The pipeline combines light sheet-based fluorescence microscopy of optically cleared spheroids with automated nuclei segmentation (F score: 0.88) and concepts from graph analysis and computational topology. Incorporating cell graphs and alpha shapes provided more than 30 features of individual nuclei, the cellular neighborhood and the spheroid morphology. The application of our pipeline to a set of breast carcinoma spheroids revealed two concentric layers of different cell density for more than 30,000 cells. The thickness of the outer cell layer depends on a spheroid’s size and varies between 50% and 75% of its radius. In differently-sized spheroids, we detected patches of different cell densities ranging from 5 × 10⁵ to 1 × 106 cells/mm³. Since cell density affects cell behavior in tissues, structural heterogeneities need to be incorporated into existing models. Our image analysis pipeline provides a multiscale approach to obtain the relevant data for a system-level understanding of tissue architecture.

Three-dimensional in vitro cell cultures more closely resemble the cellular microenvironment of cells in tissues than two-dimensional monolayer cultures¹. Compared to real tissues, they excel with well-defined experimental conditions. Even simple model systems such as monotypic spheroids² or organoids³ that show a moderate complexity, provide an adequate and reproducible characterization. Spheroids are three-dimensional multicellular clusters that form through cell aggregation and cell proliferation. With diameters of more than 400–500 μm, they develop a concentric cell layering, in which a necrotic core is surrounded by a layer of quiescent cells and an outer rim of proliferating cells⁴. Many spheroids display properties characteristic of their ancestral tissue such as beating cardiomyocyte spheroids⁵ or aggregates of mouse embryonic stem cells that exhibit axis elongation⁶. Due to their high potential, the applications of spheroids range from fundamental questions underlying cell differentiation and cancer biology to drug discovery and drug response studies⁷.

All these applications depend on the properties of individual cells in a spheroid and all means to retrieve the properties rely on spheroid disintegration or the use of rather small spheroids of less than 200 μm in diameter, which lack the prominent concentric layering and central necrosis. However, morphometric measurements in intact, differently-sized spheroids are needed⁸. Based on histological sections of spheroids, Jagiella et al. demonstrated the high potential of retrieving insight into spheroid growth from image-based modelling⁹. Automated image analysis pipelines are required to extend such an analysis to multiple cell lines. Mathematical models of spheroids have shown that changes in the state of each individual cell¹⁰ must be monitored when studying spheroids. This requires the quantification of the three-dimensional cell environment¹¹, since cell properties are affected by cell density¹², contact inhibition or mechanical pressure⁹¹³¹⁴¹⁵.

Phase contrast imaging provides a global but two-dimensional picture of spheroids¹⁶. Histological sectioning of spheroids allows the quantification of individual cells in a single two-dimensional plane⁹. Confocal and light sheet-based fluorescence microscopy (LSFM) allow the imaging of intact three-dimensional spheroids at the level of individual cells¹⁷¹⁸. For large spheroids, optical clearing methods increase object transparency by achieving refractive uniformity and increase the penetration depth¹⁹. The combination of optical clearing and LSFM has been used to analyze drug induced cell death in tumor spheroid core regions²⁰. Three-dimensional microscopy provides image stacks of fluorescently labelled cell nuclei valuable for the quantitative characterization of the three-dimensional environment in spheroids²¹. Furthermore, nuclei properties correlate with cell proliferation activity²² and have shown diagnostic value in oncology²³.

A large pool of powerful cell nuclei segmentation methods exists, including iterative thresholding²⁴, level sets²⁵, graph cut²⁶, gradient flow tracking²⁷, lines-of-sight²⁸ or watershed methods²⁹. Identifying the approach that is appropriate for a wide variety of datasets, fast and robust with respect to high cell density, as well as variations in cell nuclei volume, shape and dye distribution, has become a major challenge in image analysis. A system-level quantification requires the integration of image segmentation with powerful concepts from other fields such as graph theory, computational topology and spatial statistics. For example, the two-dimensional spatial cell network in breast tissue³⁰, malignant glioma³¹ or Hodgkin’s lymphoma³² has been characterized by cell graphs. For the implementation of complex image analysis pipelines, software environments such as Mathematica (Wolfram Research Inc.) or Matlab (MathWorks Inc.) offer comprehensive platforms that integrate well-established image analysis algorithms with a variety of techniques from other computational fields such as graph theory, statistics and computational topology. These platforms can be further extended by integrating packages such as the Insight Segmentation and Registration Toolkit (ITK)³³, the Visualization Toolkit (VTK)³⁴, Fiji³⁵ and R³⁶.

We developed a robust, multiscale approach for the characterization of large spheroids. Our approach includes three-dimensional cell culture, optical clearing, LSFM imaging and system-level image analysis. Algorithms from graph theory and computational topology complete the segmentation of cell nuclei. The integration of the Laplacian of Gaussian filter into a marker-controlled watershed algorithm provides a robust and accurate cell nuclei segmentation with an F score of 0.88. As a reference, our previous detailed analysis of available tools yielded F scores of at most 0.8²⁸. We extended cell graphs to analyze the three-dimensional spatial cell network and introduced the alpha shape as a geometrical model of spheroids. The image analysis pipeline was implemented in Mathematica and a user interface is provided.

We applied our image analysis pipeline to characterize size-dependent differences in the internal morphology of spheroids generated from breast cancer cells. Our results revealed the heterogeneity of three-dimensional superstructures that could not have been investigated so far. We detected the concentric cell layering for total cell numbers above 30,000 cells. The relative thickness of the outer region decreases from 75% to 50% of the spheroid radius with increasing cell number. The cell density in spheroids varies between 5 × 10⁵ and 1 × 10⁶ cells/mm³. Our image analysis pipeline provides the first quantitative representation of the three-dimensional cell environment in intact, differently-sized spheroids.

Results

The combination of optical clearing and LSFM provides insight into the structure of large multicellular spheroids

We applied the complete pipeline to a set of sixteen T47D spheroids that were seeded from 500 to 10,000 cells, developed for two weeks, optically cleared and, finally, imaged in toto with LSFM³⁷. This resulted in one image stack per dataset with a homogenous signal to noise ratio throughout the entire specimen (Fig. 1). Spheroid diameters range from 150 μm to more than 500 μm.

Robust and accurate segmentation of cell nuclei in large, densely packed multicellular spheroids

An accurate cell nuclei segmentation is essential to characterize spheroids quantitatively. However, imaging large samples in toto limits the achievable lateral and axial resolution. Furthermore, the high variability of cell nuclei volume and staining in spheroids impedes accurate cell nuclei segmentation. In combination with high cell density, this results in apparently touching cell nuclei in the images that are difficult to separate (Fig. 1, magnification). Our image analysis pipeline identifies cell nuclei accurately and extracts the features required for the characterization of spheroids (Fig. 2).

We use local adaptive thresholding as an initial segmentation of candidate regions. An adaptable three-dimensional, multiscale Laplacian of Gaussian (LoG) filter detects marker points reliably and reproducibly (Supplementary Figure 1). The maximal magnitude of the LoG response is achieved when the scale of the LoG matches the size of the cell nucleus. By incorporating multiple scales, the marker point detection algorithm is more robust towards variations in cell nuclei volume. We measured minimal and maximal radii rmin and rmax of cell nuclei in the datasets (three and six voxels, respectively) to determine the scales σmin and σmax of the LoG filter (see Methods section). The resulting marker points are used for a marker-controlled three-dimensional watershed segmentation of the cell nuclei (Supplementary Figure 2). The number of marker points determines the number of objects that are extracted by the watershed algorithm. The processing time for cell nuclei segmentation was less than one hour per dataset.

The segmentation performance was evaluated for three different regions I, II and III (Fig. 3a). The regions were chosen to represent the variability in cell nuclei morphology and cellular density (Fig. 3b, second column). In region I, the cell density is high, cell nuclei are morphologically very diverse and apparently touch each other. Region II contains patches of high and low cell density and exhibits a high morphological diversity. Region III contains predominantly small and spherical cell nuclei at a low density. We first evaluated the steps of the segmentation for the three regions qualitatively. The initial segmentation separated image regions that contain cell nuclei from the background. However, this step was not capable of separating clusters of apparently touching cell nuclei (Fig. 3b, third column, yellow arrows).

The marker point detection accurately identified the locations of cell nuclei and is robust in terms of size, shape and intensity variations (Fig. 3b, fourth column). Apparently touching cell nuclei were separated using the marker-controlled watershed algorithm, which was initialized at marker points identified in the marker point detection step. Borders between apparently touching cell nuclei were reconstructed and resulting cell nuclei shapes are now adequately represented (Fig. 3b, fifth column). In summary, cell nuclei in all three tested regions were accurately identified in the final segmentation. To confirm our findings of the visual inspection, we conducted a quantitative performance evaluation of the cell nuclei segmentation (Supplementary Table 1). We generated ground truth datasets (GT) for the three regions (shown in Fig. 3b) by visually extracting the locations of all cell nuclei from dataset L3 (Supplementary Table 4) with a custom program. To assess the segmentation performance, the number of correctly detected (true positives), falsely detected (false positives) and undetected cell nuclei (false negatives) were determined by matching the centroids in the GT with those resulting from the segmentation. We then derived the metrics recall, precision and F score. Out of 252 cell nuclei in the GT in region III, 216 cell nuclei were correctly identified, 36 were not detected and 9 were falsely detected (recall: 0.86, precision: 0.96, F score: 0.91). For region II, 230 out of 269 were correctly identified, 39 were not detected and 32 were falsely detected (recall: 0.86, precision: 0.88, F score: 0.87). In region I, the segmentation performance was 198 true positives out of 233 cell nuclei, 35 false negatives and 24 false positives (recall: 0.85, precision: 0.89, F score: 0.87). Overall, the number of cell nuclei detected by the segmentation is similar to that in the GT in all three regions. In total 110 false negatives and 65 false positives were obtained in all regions, indicating high precision and a low rate of under-segmentation (average recall: 0.86, average precision: 0.91, average F score: 0.88). The increased number of false positive detections in regions II and III indicates that cell nuclei segmentation in these regions is rather difficult. In region III, where cell nuclei have a similar shape and are less densely packed, the performance was slightly better. The obtained F scores indicate high efficiency and robustness of the proposed cell nuclei segmentation. Our previous, detailed performance evaluation of established nuclei segmentation approaches for T47D spheroids²⁸ yielded F scores less or equal to 0.8.

The range parameters for the initial segmentation and marker point detection are the main parameters that need to be adjusted for cell nuclei segmentation. The average diameter of the cell nuclei provides a good starting value for these parameters. A user interface allows easy fine-tuning of the parameter values (Supplementary Figures 3–5). For the underlying datasets, the same set of empirically determined parameter values was used (Supplementary Table 2, Supplementary Figure 5). However, we observed a robustness of the segmentation results with respect to changes in the parameter values.

Quantitative features characterizing spheroids

We distinguish low level features of the cell nuclei that are directly obtained from the segmentation and higher order features of the cell neighborhood and the whole spheroid. The output of our image analysis pipeline is a comprehensive set of features (Supplementary Table 3). The alpha shape was computed for the cell nuclei centroids and the boundary region was extracted as the spheroid surface (Fig. 4a). A value of 90 voxels for alpha led to a smooth approximation of the surface (Supplementary Figures 6 and 7). From the alpha shape and the extracted surface, we determined the volume, surface area, centroid and the minimal distance of the centroid to the surface.

We observed a strong positive correlation between the number of seeded cells and the spheroid volume (Pearson’s correlation coefficient: 0.95). In addition, the spheroid volume was found to be proportional to the final number of cells and we obtained an increase in volume of 1,127 μm³ per cell (Fig. 4a). This indicates a robust linear relation that exists despite the differences in the internal structure of small, medium and large spheroids. Based on automated cluster analysis (partitioning around medoids with squared Euclidean distance) of spheroid volume and cell number, we separated the datasets into three groups of small (n = 9), medium-sized (n = 3) and large (n = 4) spheroids. We determined the mean and standard deviation of the cell numbers (small: 4,977 ± 2,766, medium: 25,582 ± 454, large: 34,742 ± 2,941) and the mean and standard deviation of the spheroid volume (small: 6.0 × 10⁶ ± 3.1 × 10⁶ μm³, medium: 2.7 × 10⁷ ± 0.04 × 10⁷ μm³, large: 4.0 × 10⁷ ± 0.2 × 10⁷ μm³, Supplementary Table 4).

The volumes of the identified cell nuclei were not normally distributed (Fig. 4b). All three distributions exhibit a peak around 250 μm³ and an additional shoulder at around 450 μm³. Focusing on cell nuclei with volumes between 300 and 600 μm³, we did not observe a pattern in their spatial localization (Supplementary Figure 8). All distributions were asymmetric with positive skewness (small: 1.77, medium: 1.71, large: 1.71) and high kurtosis values (small: 8.48, medium: 7.63, large: 7.70). We obtained the median and the median absolute deviation (MAD) values for each spheroid group and found that the median cell nucleus volume in medium and large spheroids was slightly lower compared to small spheroids (small: 267 ± 68 μm³, medium: 234 ± 51 μm³, large: 225 ± 61 μm³). This information does not provide enough detail to analyze differences between the spheroid groups.

Thus, we computed the distance of each cell nucleus to the surface of the spheroid to determine its relative position in the spheroid. For a comparison between different datasets, the location is defined as the normalized shortest distance to the surface (NDS) such that cell nuclei with an NDS of 0 are at the surface and those with an NDS of 1 are in the center of a spheroid. We extended cell graphs to three spatial dimensions and computed two graph representations that capture the spatial arrangement of cells within spheroids. In these graphs, vertices represent the cell nuclei and pairs of vertices are connected by edges representing a pairwise neighborhood relationship. We distinguish between the proximity cell graph (PCG), in which edges are solely created according to the Euclidean distance and the Delaunay cell graph (DCG), in which edges are furthermore only possible between vertices that are connected by a line in the corresponding Delaunay triangulation. For both graphs, we used a distance threshold of 40 voxels for generating the edges. We employed these two graphs to analyze the internal morphology of spheroids. We determined the number of neighbors (i.e., the vertex degree) and the distance to neighbors (i.e., the weights of all incident edges) for each vertex in the graphs as cell density measures. The number of neighbors in the PCG can be interpreted as a measure of local cell density because all cells within a certain distance share the neighborhood relationship. We refer to the number of neighbors for a vertex in the PCG as the local cell density in cells/unit volume (cells/u.v.), where a unit volume corresponds to 65450 μm³.

Local cell density features reveal structural heterogeneity at multiple scales in spheroids

We inspected all features of the cell nuclei as a function of depth in the spheroid to check for variations along the radial direction. Low level cell nuclei features directly extracted from the segmentation did not show any variation along the radial direction of the spheroids. However, higher order features derived from the cell graphs that characterize local cell density varied along the radial direction in medium and large spheroids (Fig. 5a, Supplementary Figure 9). A prerequisite for any statistical analysis of cell aggregates is that measured features show a significant deviation from randomness³²³⁸. To assess the randomness of cellular arrangement, we placed cells randomly into the alpha shape of each dataset, computed the cell density features and compared them with those found for the real datasets. Thus, we plotted the mean cell density as a function of depth and compared it to the random cell position (RCP) model³⁸ (Fig. 5b). We found that the curves for medium and large spheroids strongly deviate from the RCP model, whereas in small spheroids a similar curve is obtained. Thus, the spatial distribution of cells in medium and large spheroids deviates from randomness, whereas the spatial distribution in small spheroids could be generated by randomly positioning the cells. Consequently, we further analyzed the features capturing local cell density in medium and large spheroids. Cells at the surface have only few neighboring cells (mean cell density ≈ 45 cells/u.v.). At 0.5 NDS, the cell density decreases from 65–40 cells/u.v. in large spheroids, whereas a similar decrease from 60–45 cells/u.v. is observed at 0.75 NDS in medium spheroids (Fig. 5a).

To compare our results to the known concentric cell layering in larger spheroids⁴, we manually subdivided medium and large spheroids into three distinct regions surface, outer and core. For the separation of surface and outer region we take the transition point between the initial rise and the approximately constant region, resulting in a threshold of 0.1 NDS for medium and large spheroids. For the separation of outer and core region we take the transition point between the approximately constant region and the decreasing part resulting in 0.75 NDS for medium and 0.5 NDS for large spheroids. For the mean cell density, differences were detected for outer and core regions in both spheroid groups (Fig. 5c). In medium spheroids, the mean cell density between outer and core region differs by 9.4% (outer: 61 cells/u.v., core: 56 cells/u.v.), whereas in large spheroids, the difference is 17.7% (outer: 65 cells/u.v., core: 53 cells/u.v.). To cross-check these findings, we colored the segmented cell nuclei according to the cell density (Supplementary Figure 10). Consistent with the quantification, low cell density was observed in the central and surface regions of the spheroids, while the cell density is higher for cells between these two regions (Fig. 5d). Within the identified regions, the cell density was not homogeneous and we observed patches of high cell density. This result also extends to small spheroids (Supplementary Figure 10). To evaluate the contribution of cell divisions to cell density, we manually identified occurring cell divisions in three representative datasets (S9, M2 and L3). We found between one (S9) and five (L3) cells that were at the end of the mitotic phase. For medium and large spheroids, we could not identify any correlation between the measured cell density and the nuclei volume (average Pearson’s correlation coefficient −0.23) or shape (average Pearson’s correlation coefficient −0.11). Similar to the results obtained from the PCG, the mean distance to neighbors determined from the DCG was constant up to 0.75 NDS in medium spheroids and increased to 15.5 μm in the core. In large spheroids, a similar increase was already observed at 0.5 NDS (Fig. 5e). Based on the quantitative measures, we determined the average radius of the spheroids, the average thickness of the outer region and the volume of the outer region (Fig. 5f, Supplementary Figure 11). On average, medium spheroids have a radius of 191.3 μm and large spheroids have a radius of 213.7 μm. As indicated by the previous results, the thickness of the outer region in medium spheroids is larger than in large spheroids (143.5 μm in medium, compared to 106.8 μm in large spheroids) and the volume of the outer region in medium spheroids is smaller than in large spheroids (0.028 mm³ in medium and 0.033 mm³ in large spheroids, Supplementary Figure 11). We again cross-checked by highlighting the detected border between outer and core region in the raw data (Fig. 5g). We found that the quantitatively identified location of the border fits with the visually observable boundary.

Discussion

Three-dimensional cell aggregates continue to become more and more important model systems in basic cell as well as in translational biology. Common experiments in three-dimensional cell biology rely on chemical³⁹ or mechanical⁴⁰ perturbations of spheroids. Examples include measuring the influence of mechanical confinement on the internal morphology of spheroids or the influence of chemical compounds on cell differentiation, proliferation and viability. These approaches are completed by our method, which enables the quantitative characterization of perturbation effects on spheroids.

Despite recent advances in three-dimensional cell cultures and fluorescence microscopy, any quantitative analysis of spheroids remains complex. We developed a pipeline including three-dimensional cell culture, optical clearing, LSFM imaging and three-dimensional image analysis. Apart from a detailed characterization of spheroids, our image analysis also results in a reduction of large images to simple lists. The number of extracted nuclei features goes well beyond those available in most existing programs²¹⁴¹⁴²⁴³⁴⁴. The output comprises measurements obtained from the segmentation extended by features obtained from graph and topological analysis.

The combination of optical clearing and LSFM resulted in full penetration depth, homogeneous intensity distribution and good axial resolution. We obtained the first comprehensive dataset of high quality images for subsequent quantitative characterization of spheroids at the single cell level. Our image analysis pipeline completes cell nuclei segmentation by powerful methods from graph theory and computational topology.

The cell nuclei segmentation is objective, robust and exceeds the performance of well-established programs such as FarSight and CellProfiler. It is insensitive towards intensity variations and capable of separating apparently touching cell nuclei in regions of different cellular density. The initial local segmentation is less sensitive towards heterogeneous intensity distributions than global thresholding⁴³⁴⁴⁴⁵. The multiscale LoG approach is capable of reliably detecting marker points for irregular nuclear shapes, varying intensity distributions and in regions of high cell density⁴²⁴³⁴⁵. In conjunction with results from recent studies⁴³⁴⁵, we conclude that the identified combination of methods is a promising candidate for a universally applicable cell nuclei segmentation approach.

We demonstrated that alpha shapes provide a good geometrical model of arbitrarily shaped cell aggregates, based on the positions of cells and a suitable parameter value for alpha. Compared to the commonly used spherical harmonics approach⁴⁶⁴⁷, the alpha shape excels in applicability (own unpublished results). The computational effort is low and they provide an accurate geometrical model of solid three-dimensional objects of arbitrary shape.

Graph based approaches give a formal and well-researched description of complex, multi-dimensional networks. Cell graphs inherently contain the information to compute global and local features that define the cell aggregate topology and relationships between individual cells³². To investigate the local cell structures in the spheroid, we extended cell graphs to three-dimensional spatial networks and implemented two cell graphs: a purely distance-based variant (proximity cell graph, PCG) and a variant that is based on Delaunay triangulation (Delaunay cell graph, DCG). By incorporating the alpha shape surface, these graphs capture structural patterns as a function of depth in spheroids. PCG and DCG model different modes of cell neighborhood. The PCG provides a measure of local cell density and connects all cells that interact across long ranges. Neighborhood in the DCG is an approximation of which cells are in direct contact with each other. For the PCG, the number of edges and hence the computation time increases with the distance threshold. The DCG is less sensitive to the chosen distance threshold, since the number of edges is restricted by the Delaunay triangulation. Since the Delaunay triangulation is the dual form of the Voronoi tessellation, an approximation of cellular shape can be readily obtained from the DCG⁴⁴⁴⁸. In this study, both cell graphs identified the same structural patterns in spheroids. In the future, features extracted from the segmentation can be incorporated into the cell graphs to derive higher order information of the underlying system, in particular, spatial correlation of cell morphology, cell type or expression profiles.

We applied the complete pipeline to sixteen differently-sized T47D spheroids using identical parameter values for the image analysis. Our pipeline provides quantitative data on cell location for advanced statistical analyses³⁷ and cell nuclei volume, cell number and the spheroid diameter for validation and refinement of mathematical models for spheroid growth dynamics¹³. Our approach differs from the classical approaches for studying spheroid growth dynamics. We generated spheroids from fixed numbers of seeded cells and grew them for the same period. This provides the basis for a comparison of the properties of differently-sized spheroids. We found a linear relationship between spheroid volume and cell number despite the differences in cell density distribution between small, medium and large spheroids. Hence, for a clear distinction between spheroids of different internal morphology it is insufficient to analyze such global spheroid properties. By extending the classical measures with results from cell graphs and alpha shapes, we obtained an objective method to detect local structural properties and the boundary between the outer and core region in spheroids. An extensive core region in the breast carcinoma spheroids arises for a total cell number of at least 30,000 cells. Furthermore, we revealed a structural inhomogeneity that goes well beyond the expected concentric layering. The diameter of the core region is not proportional to the diameter of the whole spheroid. In larger spheroids, the outer region is relatively thinner (50% of the spheroid radius) but has a larger volume than in medium-sized spheroids (75% of the spheroid radius). The thickness of the outer region depends on nutrient and oxygen uptake² that vary with spheroid size⁴⁹.

Within the spheroids of all sizes, the cell density varies between 35 and 70 cells/u.v. and patches of higher and lower cell density occur in the outer region. We could not identify any correlation between the cell density differences and cellular growth properties such as cell nuclei volume and the occurrence of cell divisions. The low number of cell divisions matches a cell cycle duration of 72 hours for T47D cells as measured in our laboratory (unpublished). We expect the observed heterogeneity in the cell density to have an impact on individual cells. Studies based on two-dimensional cell cultures have revealed that cell density affects cellular behavior and function⁵⁰ including cell differentiation⁵¹ and the response to compounds⁵². These findings also extend to real tissues⁵³. Hence, the heterogeneity of three-dimensional cell cultures has to be integrated into existing models.

Our approach can be readily applied to analyze spheroid growth dynamics. For example, spheroids could be grown from the same number of seeded cells for different periods to study the development of the internal spheroid morphology as a function of time. This approach could be extended to analyze growth dynamics that exceed the culturing period by incorporating the growth curves for different numbers of seeded cells into one curve. For our data, this approach indicates an exponential growth for the first 28 days, followed by a linear growth (results not shown) in agreement with previous work¹³. Furthermore, these approaches could be extended by antibody staining optimized for spheroids⁵⁴ to link the internal spheroid morphology to cellular properties (e.g. apoptotic, necrotic, proliferative). In particular, previous results have indicated a positive correlation between cell nucleus volume and cell proliferation²². We found a uniform distribution of large cell nuclei within spheroids of all sizes (Supplementary Figure 8) that interestingly deviates from the distribution of proliferating cells in recent images of histological sections stained with a proliferation marker⁹. More detailed investigations could shed light on whether this is a cell type specific phenomenon.

Our approach is not restricted to spheroids. Provided high-quality images with homogeneously stained nuclei of approximately convex shape, our image analysis can be applied to characterize any three-dimensional multicellular system at multiple scales, ranging from the single cell level to the cell microenvironment and the whole system (Supplementary Figure 12). The complete pipeline can be extended to specimens of several mm in diameter that would comprise at least several 100,000 cells. We conclude that studies in three-dimensional cell biology need to incorporate quantitative measurements at the level of the single cell, the cell neighborhood and the whole system to draw physiologically relevant conclusions (Fig. 6). A multiscale characterization can only be obtained by integrating powerful algorithms from other computational fields into traditional image analysis approaches.

Methods

Spheroid preparation and imaging

Spheroid preparation

Spheroids were formed by liquid overlay from non-invasive T47D breast cancer cells. Cell suspensions containing seed numbers of 500, 1,000, 2,000, 5,000 or 10,000 T47D breast carcinoma cells were placed in U-well plates coated with a non-adhesive surface (Thermo Scientific, Nunc HydroCell 96 Well, 174908) that facilitates spheroid formation. The cells were cultured for twelve days in order to obtain mature spheroids. Subsequently, spheroids were chemically fixed, washed, stained for cell nuclei with the fluorescent dye Draq5 and optically cleared with BABB¹⁹ to minimize light scattering.

Image acquisition

Image stacks of optical sections of entire T47D spheroids were acquired with a monolithic Digital Scanned Laser Light Sheet-based Fluorescence Microscope (mDSLM) as described previously⁵⁵. We define the illumination axis as X, the detection axis as Z and the axis orthogonal to X and Z as Y. The excitation light source is a 638 nm laser operated at 1.6 mW and focused with a Carl Zeiss Epiplan-Neofluar 2.5x/NA 0.06 objective lens. A Carl Zeiss N-Achroplan 20x/NA 0.5 objective lens, a 680/60 bandpass filter, a regular tube lens and a high-resolution Andor Neo 5.5 sCMOS camera with a pixel pitch of 6.5 μm are used in the emission path. Image stacks were acquired with an axial pitch of 1.29 μm. The cropped regions of interest result in raw image stack sizes ranging from 230 Megabytes to 2 Gigabytes.

Cell nuclei segmentation

Pre-processing

Isotropic voxels are obtained by rescaling the Z dimension and interpolating missing planes (ImageZScalingFactor). This results in isotropic voxels with a pitch of 0.325 μm. For computational efficiency, the resulting three-dimensional image is resized by a factor of 0.5 (ImageScalingFactor). In the resized image, the voxels are isotropic with a pitch of 0.65 μm. In the following, the raw three-dimensional image is termed fr(x, y, z).

Initial cell nuclei segmentation

Noise in fr(x, y, z) is reduced by the convolution with a three-dimensional Gaussian kernel of range 3 × 3 × 3 (NucleiFilterRange). The average background intensity tglobal of the convolved image is determined by Otsu’s global thresholding algorithm⁵⁶. In the initial cell nuclei segmentation step, local thresholding is applied per sectional plane along the dimensions X, Y and Z and the resulting binary images are multiplied to obtain the initial segmentation. For each pixel the local threshold tlocal is determined by

where mlocal corresponds to the mean intensity measured in a range of 25 × 25 pixels (NucleiThresholdRange) and tglobal is the determined background intensity. The factor γ is a parameter of the cell nuclei segmentation that controls the impact of the determined background intensity (NucleiBackgroundFactor). For all datasets, γ was set to 0.25. Each pixel with an intensity greater than tlocal is set to 1, all others are set to 0. Holes, i.e. regions that are falsely detected as background because of minor segmentation errors, are removed. Therefore, the initial segmentation is inverted to obtain the complement image. Holes are identified as small foreground regions and those smaller than a predefined threshold of 250 voxels are removed. This threshold corresponds to about one quarter of the mean volume of a cell nucleus and was determined prior to segmentation. The image is inverted again to obtain the hole-corrected fb(x, y, z).

Marker point detection

For the decomposition of connected nuclei in fb(x, y, z), we use a three-dimensional marker-controlled immersion watershed algorithm. The marker points are obtained by the multiscale Laplacian of Gaussian (LoG) blob detection algorithm. Thereby, a blob refers to an approximately convex region of a three-dimensional image, in which the intensities vary within a sufficiently small range of values. First, the raw image fr(x, y, z) is inverted, resulting in fi(x, y, z). The LoG corresponds to first convolving fi(x, y, z) by a Gaussian kernel g(x, y, z; σ). Thus, we get

where σ is the standard deviation (scale) of the Gaussian kernel g. Then the Laplacian operator ∇² is applied to the convolved image to obtain the LoG response. In our approach, the inverted image fi(x, y, z) is processed at multiple scales . The minimal and maximal scales σmin and σmax are determined using the relationship

between radius r of a blob-like object and the scale σ of the LoG. We measured the minimal (rmin) and maximal (rmax) radius of cell nuclei a priori in the images (NucleiSeedDetectionMinRadius and NucleiSeedDetectionMaxRadius) and used the relationship above to compute σmin and σmax. To achieve scale-invariance, the LoG response at scale σ is normalized by multiplication with σ³ such that the LoG response at scale σ is given by

For computational efficiency, we iteratively compute the maximum response over scales given by

The maximum response is further processed with a maximum transform for detection of extended maxima. The obtained binary image is multiplied with the initial segmentation fb(x, y, z) to discard extended maxima detected in the background. Extended maxima that are in close proximity to each other are merged by increasing the size of the marker points using a morphological dilation operator with a round structuring element of two voxels radius (NucleiSeedDilation). The resulting binary image fm(x, y, z) specifies the marker points that are used to initialize the subsequent watershed algorithm.

Watershed-based decomposition of cell nuclei clusters

We use the immersion watershed algorithm to achieve a decomposition of connected cell nuclei clusters. The algorithm starts at the marker points in fm(x, y, z) as sources and iteratively immerses the inverted image fi(x, y, z) according to the intensity. Watersheds are built when different sources meet during the immersion process. A partitioning of the image fi(x, y, z) into labelled components (1 to number of components) and watersheds (0 s) is obtained. However, the immersion process does not necessarily stop at the border of a cell nucleus. Thus, we multiply the resulting matrix with the initial segmentation fb(x, y, z). Objects are selected by incorporating lower and upper volume thresholds (NucleiMinCount and NucleiMaxCount). Based on the measurements of minimal and maximal radii of cell nuclei, the cell nucleus volume is approximated as a sphere with equivalent radius and we obtained an approximate lower threshold of 250 voxels and an upper threshold of 42,500 voxels. As the final result of the segmentation, we obtain the matrix w(x, y, z). Each labelled component in w(x, y, z) represents a cell nucleus found by the segmentation.

Evaluation of segmentation performance

We generated three ground truth (GT) datasets by cropping three sub-regions of 100 × 100 × 100 voxels out of the downsized raw images. The centroids of all cell nuclei within the sub-regions were visually identified with a custom program. Based on the generated GT, we determined the number of correctly detected cell nuclei (true positives, TP), the number of cell nuclei that were falsely detected by the segmentation (false positives, FP) and the number of cell nuclei in the GT that were not detected by the segmentation (false negatives, FN). To compute these numbers, we used the following algorithm: if exactly one centroid of the segmentation is found within a spherical neighborhood of twelve voxels of a centroid in the GT, we count it as TP and delete it from the list. If more than one centroid is found within this neighborhood range, the closest one is considered as TP. Based on the number of TP, the FP and FN were obtained using FP = NSC−TP and FN = NGT − TP, where NSC is the number of centroids determined by the segmentation and NGT is the number of centroids in the GT. Based on these measurements we derived the metrics recall, precision and F score²⁸ with values ranging from 0 (worst performance) to 1 (optimal performance):

Feature computation and randomness of the spatial cell distribution

Basic cell nuclei features

Based on the segmentation, we extract intensity-related and morphological features of each identified cell nucleus. The morphological features include the volume in number of voxels, centroid, intensity-weighted centroid, mean, minimum and maximum distance to the centroid, number of voxels on the surface and the bounding box. We use principal component analysis to determine the principal directions and the extensions along these directions as a measure of shape and orientation of the cell nucleus. Intensity-related features include the mean, minimum and maximum intensity and the intensity standard deviation. All intensity-related features are rescaled to the interval [0, 1].

Surface approximation

We approximate the surface of the spheroid using the alpha shapes approach based on the set P of cell nuclei centroids. To discriminate against false positive detections, all points in P that are within a distance of 20 voxels are connected (OutlierDistanceThreshold). The largest connected component is regarded as the cell aggregate and remaining components are removed. The surface of the cell aggregate is then approximated using Edelsbrunner’s algorithm for alpha shapes⁵⁷ with parameter alpha (Alpha). The alpha shape approach is based on a Delaunay triangulation and works well for point sets P of high density and uniform distribution. If the corresponding parameter alpha approaches 0, the alpha shape converges to P; whereas if alpha approaches ∞, the alpha shape converges to the convex hull of P. We set alpha to 90 voxels for all datasets. The obtained alpha shape is used to compute the volume, extract the surface and compute the surface area of the spheroid. Having obtained the surface of the spheroid, the relative position of each cell nucleus is determined by computing its distance to the surface.

Cell graphs

We derive two graph representations from the cell nuclei centroids representing the spatial distribution of cells within the cell aggregates. The proximity cell graph is given by PCG(V, EPCG) where V is the vertex set (i.e. the cell nuclei centroids) and EPCG is the edge set of the graph. In the PCG, cells are neighbors if they are closer than a certain distance (EdgeDistanceThreshold). Thus, we obtained an edge (u, w) between two vertices u and w if the Euclidean distance between u and w is less than a predefined threshold of 40 voxels. The Delaunay cell graph is given by DCG(V, EDCG) where V is the vertex set and EDCG is the edge set of the graph. The DCG graph is constructed based on a Delaunay triangulation to approximate which cells are in physical contact. An edge (u, w) was created between two vertices u and w if the corresponding points are connected by a line in the Delaunay triangulation and the Euclidean distance between u and w was less than a predefined threshold (EdgeDistanceThreshold) of 40 voxels. Consequently, the edge set EDCG is a subset of EPCG and the vertex set V for both graphs is identical. Edge weights in both graphs are equal to the Euclidean distance between the corresponding vertices. We extract the degree of each vertex as the number of neighbors. Further, the minimum, maximum, mean and standard deviation of the distance to neighbors for each vertex v is given by the weights of all edges incident to v.

Random cell position model

We compared the generated cell graphs to those of a mathematical random cell position (RCP) model. In this model, we made the following assumptions: (1) cells are distributed uniformly within a spheroid. (2) Cell nuclei are represented as non-overlapping spheres with positions drawn from a uniform random number distribution. The radius of the spheres was set to the median cell nucleus radius (~6 voxels) from all datasets. (3) Cell nuclei positions are restricted to the volume of the spheroid (i.e. the alpha shape) and are not allowed outside this volume. (4) The number of randomly generated cell nuclei is the same as the number of cell nuclei determined for the real spheroid. For each dataset, we performed ten Monte Carlo simulations of the RCP model and obtained the proximity and Delaunay cell graphs in the same way as for the real datasets. In other words, we generated an exact analogue of each cell aggregate.

Implementation details

Implementation

The image analysis pipeline was developed in Mathematica version 10.2 (Wolfram Research, Inc.). A package containing the custom code of the image analysis pipeline, a user interface and two example datasets are available for download at http://www.physikalischebiologie.de/downloads. Computations were conducted on a workstation that comprises two six-core CPUs (X5650, Intel Corporation), 96 Gigabyte DDR3 memory, running Windows Server 2012 R2. Cell nuclei segmentation and post-processing took between five minutes and less than one hour per stack. The current implementation is parallelized at the initial cell nuclei segmentation step. For all subsequent steps, the complete image needs to fit into the main memory.

Three-dimensional reconstruction

Three-dimensional reconstructions were generated using the Arivis software (http://www.arivis.com).

Additional Information

How to cite this article: Schmitz, A. et al. Multiscale image analysis reveals structural heterogeneity of the cell microenvironment in homotypic spheroids. Sci. Rep. 7, 43693; doi: 10.1038/srep43693 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figure 1View larger versionFigure 1. Image quality of three-dimensional datasets.Three-dimensional volume rendering (first column), single planes along X-Y (second column), single planes along Z-Y (third column) and magnification (fourth column) of two spheroids of 500 (upper row, dataset S9) and 10,000 (lower row, dataset L3) seeded cells. For a complete list of datasets see Supplementary Table 4. Renderings in the first column were clipped at about the center of the spheroids and single planes were taken at the same position. Yellow boxes indicate the parts of the images magnified in the fourth column. Microscope: mDSLM. Excitation lens: CZ 5x/NA 0.16. Emission lens: CZ 20x/NA 0.50. Scale bars: 50 μm for the single planes in the second and third column, 10 μm for magnified images in the fourth column.

Figure 2View larger versionFigure 2. Illustration of the main steps of the automated image analysis pipeline.(a) The raw data comprises a stack of two-dimensional optical section images, which in this case consists of 218 planes. (b) The image stack is interpolated to obtain isotropic voxels and scaled down by a factor of 0.5. (c) Marker positions (indicated as white crosses) are identified by multiscale Laplacian of Gaussian (LoG) filtering of the raw image. (d) Regions that contain cell nuclei are separated from the image background by local thresholding of the raw image. (e) Identified marker points, the preprocessed raw image and the initial segmentation are subjected to a three-dimensional marker-controlled watershed algorithm, which separates clusters of apparently touching cell nuclei. Different colors indicate individual cell nuclei. (f) An alpha shape is constructed from the cell nuclei centroids and the boundary region is extracted as the surface (light blue). (g) Cell graphs are generated from vertices, which correspond to the cell nuclei (white spheres), and edges, which indicate the neighborhood relation between two cell nuclei (white lines). (h) Features of each individual cell nucleus, the cell neighborhood and the spheroid are extracted and stored in a tabular format.

Figure 3View larger versionFigure 3. Quality of cell nuclei segmentation in different regions of a large T47D spheroid.(a) Three-dimensional rendering of cell nuclei clipped at about the center of an exemplary spheroid (dataset L3, see Supplementary Table 4) that comprises 38,783 cells. Colored boxes indicate three regions that exhibit visually distinguishable properties. Red box: Cell nuclei in region I appear tightly packed. Green box: cell nuclei in region II exhibit diverse morphologies and the intercellular distances differ. Blue box: cell nuclei in region III of the spheroid are small, spherical and appear well separated. (b) Exemplary cell nuclei segmentation results for the three regions I, II and III. First column: sub-regions of 100 × 100 × 100 voxels were copied out of the pre-processed raw image. Second column: plane 58 of each sub-region in X-Y view. Third column: the result of the initial segmentation for plane 58. The initial segmentation accurately identifies the foreground in the image. Locations marked with yellow arrows show cell nuclei clusters that are not yet separated. Fourth column: marker points detected by the multiscale LoG filter are overlaid as yellow crosses. Note that for illustration purposes, the marker detection was performed in two dimensions. Fifth column: overlay of the final cell nuclei segmentation of plane 58 after three-dimensional marker-controlled watershed. This step effectively separates apparently touching cell nuclei (yellow arrows). Sixth column: the final segmentation result after three-dimensional marker-controlled watershed of the subregions shown in the first column. Different colors represent individual cell nuclei. Scale bar: 10 μm.

Figure 4View larger versionFigure 4. Clustering of datasets results in three groups.(a) Plot of spheroid volume versus number of cells detected in the spheroid for all datasets. The measurements are well fitted by a linear model with slope 1,118 μm³/cell, Automated clustering of the datasets according to cell number and spheroid volume, returned three groups: small (nine datasets, green ellipse), medium (three datasets, black ellipse) and large spheroids (four datasets, light blue ellipse). (b) Smoothed histogram of the cell nuclei volume distribution of small, medium and large spheroids.

Figure 5View larger versionFigure 5. Medium and large spheroids show differences in internal structure.Plot of mean cell density versus the normalized distance to the surface for small, medium and large spheroids (a) and the corresponding random cell position (RCP) models (b). The RCP models represent an exact analogue of each spheroid with the only difference that cell nuclei are randomly positioned (details are in the Materials and Methods section). The shaded regions indicate the standard error of the mean (SE). Note that in some cases, due to the small error, the shaded region is not visible. Based on the normalized distance to the surface, cell nuclei in the range [0.1, 0.75] for medium spheroids, and [0.1, 0.5] for large spheroids are assigned to the outer region, whereas cell nuclei in the range [0.75, 0.1] for medium and [0.5, 0.1] for large spheroids are assigned to the core region. (c) Mean cell density in the outer and core region for medium and large spheroids. (d) Three-dimensional rendering of segmented cell nuclei colored according to their corresponding cell density for a medium (dataset M3) and a large spheroid (dataset L2), ranging from blue (35 cells/u.v.) to red (70 cells/u.v.). For a complete list of datasets see Supplementary Table 4. Renderings were clipped at the center of the spheroids. (e) Plot of the mean distance to neighbors versus the normalized distance to the surface. The shaded regions indicate the standard error of the mean (SE). Note that in some cases, due to the small error, the shaded region is not visible. (f) Radius of medium and large spheroids and thickness of the outer region. (g) Raw single planes at the center of the spheroids shown in d with orange dashed lines indicating the boundaries of the outer region. Scale bar: 50 μm.

Figure 6View larger versionFigure 6. Three-dimensional cell biology requires a multiscale image analysis approach.The established methods for advanced three-dimensional microscopy and image segmentation need to be extended by concepts from other fields including graph theory and computational topology. This will provide a system-level understanding of tissue architecture and function.

Supplementary Material

Open In Web Browser

Acknowledgments

We thank I. Koch and J. Ackermann for critical comments and discussions on graph theoretical approaches and I. Smyrek for proof reading the manuscript. We acknowledge the funding by the Hessian LOEWE Ub-Net (speaker Ivan Dikic) and by the Cluster of Excellence Frankfurt for Macromolecular Complexes (CEF-MC, EXC 115), partially granted to the Buchmann Institute for Molecular Life Sciences (BMLS) at the Goethe Universität Frankfurt am Main by the Deutsche Forschungsgemeinschaft (DFG).

Footnotes

The authors declare no competing financial interests.

Author Contributions A.S. and S.C.F. developed the algorithms for segmentation and data analysis and wrote the manuscript. C.M. prepared the spheroids and acquired the raw image stacks. F.P. supervised the work of C.M. and contributed to the clearing and imaging experiments. EHKS supervised the work and critically revised the manuscript. All authors read and approved the final paper.

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group

PMC Copyright Notice

The articles available from the PMC site are protected by copyright, even though access is free. Copyright is held by the respective authors or publishers who provide these articles to PMC. Users of PMC are responsible for complying with the terms and conditions defined by the copyright holder.

Users should assume that standard copyright protection applies to articles in PMC, unless an article contains an explicit license statement that gives a user additional reuse or redistribution rights. PMC does not allow automated/bulk downloading of articles that have standard copyright protection.

See the copyright notice on the PMC site, https://www.ncbi.nlm.nih.gov/pmc/about/copyright/, for further details and specific exceptions.

References

Pampaloni F., Reynaud E. G. & Stelzer E. H. K. The third dimension bridges the gap between cell culture and live tissue. Nat. Rev. Mol. Cell Biol. 8, 839–845 (2007).17684528

Sutherland R. M. Cell and Environment Interactions in Tumor Microregions: The Multicell Spheroid Model. Science 240, 177–184 (1988).2451290

Marx V. Tissue engineering: Organs from the lab. Nature 522, 373–377 (2015).26085275

Kunz-Schughart L. A., Freyer J. P., Hofstaedter F. & Ebner R. The use of 3-D cultures for high-throughput screening: the multicellular spheroid model. J. Biomol. Screen. 9, 273–285 (2004).15191644

Bartholomä P. . Three-dimensional in vitro reaggregates of embryonic cardiomyocytes: a potential model system for monitoring effects of bioactive agents. J. Biomol. Screen. 10, 814–22 (2005).16234345

van den Brink S. C. . Symmetry breaking, germ layer specification and axial organisation in aggregates of mouse embryonic stem cells. Development 141, 4231–42 (2014).25371360

Ravi M., Paramesh V., Kaviya S. R., Anuradha E. & Solomon F. D. P. 3D cell culture systems: advantages and applications. J. Cell. Physiol. 230, 16–26 (2015).24912145

Bredel-Geissler A., Karbach U., Walenta S., Vollrath L. & Mueller-Klieser W. Proliferation-associated oxygen consumption and morphology of tumor cells in monolayer and spheroid culture. J. Cell. Physiol. 153, 44–52 (1992).1522135

Jagiella N., Müller B., Müller M., Vignon-Clementel I. E. & Drasdo D. Inferring Growth Control Mechanisms in Growing Multi-cellular Spheroids of NSCLC Cells from Spatial-Temporal Image Data. PLoS Comput. Biol. 12, e1004412 (2016).26866479

Powathil G. G., Swat M. & Chaplain M. A. J. Systems oncology: Towards patient-specific treatment regimes informed by multiscale mathematical modelling. Semin. Cancer Biol. 30, 13–20 (2014).24607841

Kempf H., Hatzikirou H., Bleicher M. & Meyer-Hermann M. In Silico Analysis of Cell Cycle Synchronisation Effects in Radiotherapy of Tumour Spheroids. PLoS Comput. Biol. 9, e1003295 (2013).24244120

Böttger K. . An Emerging Allee Effect Is Critical for Tumor Initiation and Persistence. PLoS Comput. Biol. 11, e1004366 (2015).26335202

Drasdo D. & Höhme S. A single-cell-based model of tumor growth in vitro: monolayers and spheroids. Phys. Biol. 2, 133–147 (2005).16224119

Montel F. . Stress Clamp Experiments on Multicellular Tumor Spheroids. Phys. Rev. Lett. 107, 188102 (2011).22107677

Cheng G., Tse J., Jain R. K. & Munn L. L. Micro-environmental mechanical stress controls tumor spheroid size and morphology by suppressing proliferation and inducing apoptosis in cancer cells. PLoS One 4, e4632 (2009).19247489

Ansari N., Müller S., Stelzer E. H. K. & Pampaloni F. Quantitative 3D Cell-Based Assay Performed with Cellular Spheroids and Fluorescence Microscopy. Methods Cell Biol. 113, 295–309 (2013).23317907

Verveer P. J. . High-resolution three-dimensional imaging of large specimens with light sheet-based microscopy. Nat. Methods 4, 311–313 (2007).17339847

Ansari N. . Quantifying the autophagy-triggering effects of drugs in cell spheroids with live fluorescence microscopy. Methods Mol. Biol. 1165, 19–29 (2014).24839016

Dodt H.-U. . Ultramicroscopy: three-dimensional visualization of neuronal networks in the whole mouse brain. Nat. Methods 4, 331–336 (2007).17384643

Wenzel C. . 3D high-content screening for the identification of compounds that target cells in dormant tumor spheroid regions. Exp. Cell Res. 323, 131–143 (2014).24480576

Bilgin C. C., Kim S., Leung E., Chang H. & Parvin B. Integrated profiling of three dimensional cell culture models and 3D microscopy. Bioinformatics 29, 3087–3093 (2013).24045773

Smolle J., Hofmann-Wellenhof R., Soyer H. P., Stettner H. & Kerl H. Nuclear size and shape parameters correlate with proliferative activity in cutaneous melanocytic tumors. J. Invest. Dermatol. 93, 178–82 (1989).2746003

Saito A. . A novel method for morphological pleomorphism and heterogeneity quantitative measurement: Named cell feature level co-occurrence matrix. J. Pathol. Inform. 7, 36 (2016).27688927

Keller P. J., Schmidt A. D., Wittbrodt J. & Stelzer E. H. K. Reconstruction of zebrafish early embryonic development by scanned light sheet microscopy. Science 322, 1065–1069 (2008).18845710

Bergeest J.-P. & Rohr K. Efficient globally optimal segmentation of cells in fluorescence microscopy images using level sets and convex energy functionals. Med. Image Anal. 16, 1436–1444 (2012).22795525

Qi J. Dense nuclei segmentation based on graph cut and convexity-concavity analysis. J. Microsc. 253, 42–53 (2014).24237576

Li G. . 3D cell nuclei segmentation based on gradient flow tracking. BMC Cell Biol. 8, 40 (2007).17784958

Mathew B. . Robust and automated three-dimensional segmentation of densely packed cell nuclei in different biological specimens with Lines-of-Sight decomposition. BMC Bioinformatics 16, 187 (2015).26049713

Cheng J. & Rajapakse J. C. Segmentation of clustered nuclei with shape markers and marking function. IEEE Trans. Biomed. Eng. 56, 741–748 (2009).19272880

Bilgin C., Demir C., Nagi C. & Yener B. Cell-Graph Mining for Breast Tissue Modeling and Classification. In 2007 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society 2007, 5311–5314 IEEE (2007).

Demir C., Gultekin S. H. & Yener B. Augmented cell-graphs for automated cancer diagnosis. Bioinformatics 21, ii7–12 (2005).16204128

Schäfer H. . CD30 cell graphs of Hodgkin lymphoma are not scale-free - An image analysis approach. Bioinformatics, doi: 10.1093/bioinformatics/btv542 (2015).

Yoo T. S. . Engineering and algorithm design for an image processing API: A technical report on ITK - The Insight Toolkit. In Studies in Health Technology and Informatics 85, 586–592 (2002).15458157

Schroeder W. J., Martin K. M. & Lorensen W. E. The design and implementation of an object-oriented toolkit for 3D graphics and visualization. Proc. Seventh Annu. IEEE Vis. ’96 1, 93–100 (1996).

Schindelin J. . Fiji: an open-source platform for biological-image analysis. Nat. Methods 9, 676–682 (2012).22743772

R Development Core Team. R: A language and environment for statistical computing. (2016).

Dini S. . Identifying the necrotic zone boundary in tumour spheroids with pair-correlation functions. J. R. Soc. Interface 13, 123 (2016).

Mokhtari Z., Mech F., Zehentmeier S., Hauser A. E. & Figge M. T. Quantitative image analysis of cell colocalization in murine bone marrow. Cytom. Part A 87, 503–512 (2015).

Friedrich J., Seidel C., Ebner R. & Kunz-Schughart L. A. Spheroid-based drug screen: considerations and practical approach. Nat. Protoc. 4, 309–324 (2009).19214182

Delarue M. . Compressive stress inhibits proliferation in tumor spheroids through a volume limitation. Biophys. J. 107, 1821–8 (2014).25418163

Lin G. . A multi-model approach to simultaneous segmentation and classification of heterogeneous populations of cell nuclei in 3D confocal microscope images. Cytom. Part A 71, 724–736 (2007).

Stegmaier J. . Fast segmentation of stained nuclei in terabyte-scale, time resolved 3D microscopy image stacks. PLoS One 9, e90036 (2014).24587204

Morales-Navarrete H. . A versatile pipeline for the multi-scale digital reconstruction and quantitative analysis of 3D tissue architecture. Elife 4, e11214 (2015).26673893

Friebel A. . TiQuant: Software for tissue analysis, quantification and surface reconstruction. Bioinformatics 31, 3234–3236 (2015).26040455

Chittajallu D. R. . In vivo cell-cycle profiling in xenograft tumors by quantitative intravital microscopy. Nat. Methods 12, 577–85 (2015).25867850

Khairy K., Reynaud E. & Stelzer E. Detection of deformable objects in 3D images using Markov-Chain Monte Carlo and spherical harmonics. Lect. Notes Comput. Sci. (including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics) 5242 LNCS, 1075–1082 (2008).

Khairy K., Foo J. & Howard J. Shapes of Red Blood Cells: Comparison of 3D Confocal Images with the Bilayer-Couple Model. Cell. Mol. Bioeng. 1, 173–181 (2008).

Schaller G. & Meyer-Hermann M. Multicellular tumor spheroid in an off-lattice Voronoi-Delaunay cell model. Phys. Rev. E 71, 51910 (2005).

Mueller-Klieser W. Method for the determination of oxygen consumption rates and diffusion coefficients in multicellular spheroids. Biophys. J. 46, 343–348 (1984).6487734

Mo J.-S., Park H. W. & Guan K.-L. The Hippo signaling pathway in stem cell biology and cancer. EMBO Rep. 15, 642–56 (2014).24825474

Gorman B. R. . Multi-scale imaging and informatics pipeline for in situ pluripotent stem cell analysis. PLoS One 9, e116037 (2014).25551762

Greene J. M., Levy D., Herrada S. P., Gottesman M. M. & Lavi O. Mathematical modeling reveals that changes to local cell density dynamically modulate baseline variations in cell growth and drug response. Cancer Res., doi: 10.1158/0008-5472.CAN-15-3232 (2016)

Mateus R. . Control of tissue growth by Yap relies on cell density and F-actin in zebrafish fin regeneration. Development 142, 2752–2763 (2015).26209644

Smyrek I. & Stelzer E. H. K. Quantitative three-dimensional evaluation of immunofluorescence staining for large whole mount spheroids with light sheet microscopy. Biomed. Opt. Express 8, 484–499 (2017).

Keller P. J. & Stelzer E. H. Quantitative in vivo imaging of entire embryos with Digital Scanned Laser Light Sheet Fluorescence Microscopy. Curr. Opin. Neurobiol. 18, 624–632 (2008).19375303

Otsu N. A Threshold Selection Method from Gray-Level Histograms. IEEE Trans. Syst. Man. Cybern. 9, 62–66 (1979).

Edelsbrunner H., Kirkpatrick D. & Seidel R. On the shape of a set of points in the plane. IEEE Trans. Inf. Theory 29, 551–559 (1983).

[Back]

Figure 1. Image quality of three-dimensional datasets.

Three-dimensional volume rendering (first column), single planes along X-Y (second column), single planes along Z-Y (third column) and magnification (fourth column) of two spheroids of 500 (upper row, dataset S9) and 10,000 (lower row, dataset L3) seeded cells. For a complete list of datasets see Supplementary Table 4. Renderings in the first column were clipped at about the center of the spheroids and single planes were taken at the same position. Yellow boxes indicate the parts of the images magnified in the fourth column. Microscope: mDSLM. Excitation lens: CZ 5x/NA 0.16. Emission lens: CZ 20x/NA 0.50. Scale bars: 50 μm for the single planes in the second and third column, 10 μm for magnified images in the fourth column.

[Back]

[Back]

Figure 2. Illustration of the main steps of the automated image analysis pipeline.

(a) The raw data comprises a stack of two-dimensional optical section images, which in this case consists of 218 planes. (b) The image stack is interpolated to obtain isotropic voxels and scaled down by a factor of 0.5. (c) Marker positions (indicated as white crosses) are identified by multiscale Laplacian of Gaussian (LoG) filtering of the raw image. (d) Regions that contain cell nuclei are separated from the image background by local thresholding of the raw image. (e) Identified marker points, the preprocessed raw image and the initial segmentation are subjected to a three-dimensional marker-controlled watershed algorithm, which separates clusters of apparently touching cell nuclei. Different colors indicate individual cell nuclei. (f) An alpha shape is constructed from the cell nuclei centroids and the boundary region is extracted as the surface (light blue). (g) Cell graphs are generated from vertices, which correspond to the cell nuclei (white spheres), and edges, which indicate the neighborhood relation between two cell nuclei (white lines). (h) Features of each individual cell nucleus, the cell neighborhood and the spheroid are extracted and stored in a tabular format.

[Back]

[Back]

Figure 3. Quality of cell nuclei segmentation in different regions of a large T47D spheroid.

(a) Three-dimensional rendering of cell nuclei clipped at about the center of an exemplary spheroid (dataset L3, see Supplementary Table 4) that comprises 38,783 cells. Colored boxes indicate three regions that exhibit visually distinguishable properties. Red box: Cell nuclei in region I appear tightly packed. Green box: cell nuclei in region II exhibit diverse morphologies and the intercellular distances differ. Blue box: cell nuclei in region III of the spheroid are small, spherical and appear well separated. (b) Exemplary cell nuclei segmentation results for the three regions I, II and III. First column: sub-regions of 100 × 100 × 100 voxels were copied out of the pre-processed raw image. Second column: plane 58 of each sub-region in X-Y view. Third column: the result of the initial segmentation for plane 58. The initial segmentation accurately identifies the foreground in the image. Locations marked with yellow arrows show cell nuclei clusters that are not yet separated. Fourth column: marker points detected by the multiscale LoG filter are overlaid as yellow crosses. Note that for illustration purposes, the marker detection was performed in two dimensions. Fifth column: overlay of the final cell nuclei segmentation of plane 58 after three-dimensional marker-controlled watershed. This step effectively separates apparently touching cell nuclei (yellow arrows). Sixth column: the final segmentation result after three-dimensional marker-controlled watershed of the subregions shown in the first column. Different colors represent individual cell nuclei. Scale bar: 10 μm.

[Back]

[Back]

Figure 4. Clustering of datasets results in three groups.

(a) Plot of spheroid volume versus number of cells detected in the spheroid for all datasets. The measurements are well fitted by a linear model with slope 1,118 μm³/cell, Automated clustering of the datasets according to cell number and spheroid volume, returned three groups: small (nine datasets, green ellipse), medium (three datasets, black ellipse) and large spheroids (four datasets, light blue ellipse). (b) Smoothed histogram of the cell nuclei volume distribution of small, medium and large spheroids.

[Back]

[Back]

Figure 5. Medium and large spheroids show differences in internal structure.

Plot of mean cell density versus the normalized distance to the surface for small, medium and large spheroids (a) and the corresponding random cell position (RCP) models (b). The RCP models represent an exact analogue of each spheroid with the only difference that cell nuclei are randomly positioned (details are in the Materials and Methods section). The shaded regions indicate the standard error of the mean (SE). Note that in some cases, due to the small error, the shaded region is not visible. Based on the normalized distance to the surface, cell nuclei in the range [0.1, 0.75] for medium spheroids, and [0.1, 0.5] for large spheroids are assigned to the outer region, whereas cell nuclei in the range [0.75, 0.1] for medium and [0.5, 0.1] for large spheroids are assigned to the core region. (c) Mean cell density in the outer and core region for medium and large spheroids. (d) Three-dimensional rendering of segmented cell nuclei colored according to their corresponding cell density for a medium (dataset M3) and a large spheroid (dataset L2), ranging from blue (35 cells/u.v.) to red (70 cells/u.v.). For a complete list of datasets see Supplementary Table 4. Renderings were clipped at the center of the spheroids. (e) Plot of the mean distance to neighbors versus the normalized distance to the surface. The shaded regions indicate the standard error of the mean (SE). Note that in some cases, due to the small error, the shaded region is not visible. (f) Radius of medium and large spheroids and thickness of the outer region. (g) Raw single planes at the center of the spheroids shown in d with orange dashed lines indicating the boundaries of the outer region. Scale bar: 50 μm.

[Back]

[Back]

Figure 6. Three-dimensional cell biology requires a multiscale image analysis approach.

The established methods for advanced three-dimensional microscopy and image segmentation need to be extended by concepts from other fields including graph theory and computational topology. This will provide a system-level understanding of tissue architecture and function.

[Back]