Spatial analysis of histology in 3D: quantification and visualization of organ and tumor level tissue environment.

Histological changes in tissue are of primary importance in pathological research and diagnosis. Automated histological analysis requires ability to computationally separate pathological alterations from normal tissue. Conventional histopathological assessments are performed from individual tissue sections, leading to the loss of three-dimensional context of the tissue. Yet, the tissue context and spatial determinants are critical in several pathologies, such as in understanding growth patterns of cancer in its local environment. Here, we develop computational methods for visualization and quantitative assessment of histopathological alterations in three dimensions. First, we reconstruct the 3D representation of the whole organ from serial sectioned tissue. Then, we proceed to analyze the histological characteristics and regions of interest in 3D. As our example cases, we use whole slide images representing hematoxylin-eosin stained whole mouse prostates in a Pten+/- mouse prostate tumor model. We show that quantitative assessment of tumor sizes, shapes, and separation between spatial locations within the organ enable characterizing and grouping tumors. Further, we show that 3D visualization of tissue with computationally quantified features provides an intuitive way to observe tissue pathology. Our results underline the heterogeneity in composition and cellular organization within individual tumors. As an example, we show how prostate tumors have nuclear density gradients indicating areas of tumor growth directions and reflecting varying pressure from the surrounding tissue. The methods presented here are applicable to any tissue and different types of pathologies. This work provides a proof-of-principle for gaining a comprehensive view from histology by studying it quantitatively in 3D.

Introduction

Tissue histology is one of the main determinants in studying and diagnosing various pathologies, including cancer. Tumor cells differ from normal cells in their appearance in tissue, and tumor cell characteristics are often utilized in cancer grading. Solid tumors may also change the structure and appearance of the surrounding tissue due to e.g. pressure, inflammation, activation of stroma, angiogenesis, and tissue stress responses they induce. Traditionally, histopathological assessments both in research and clinical settings utilize individual tissue sections in 2D, resulting in the loss of the three-dimensional context of the tissue. This has left our knowledge of the tissue context and spatial determinants in several pathologies incomplete, hampering our understanding of e.g. growth patterns of tumors in their local environment.

Currently, traditional microscopy is being increasingly substituted by digital pathology, in which tissue sections are imaged to digital high-resolution whole slide images (WSI) [20]. In addition to enabling detailed visual analysis of the histological samples, the WSIs now enable digital analysis and quantification of histological features [20], [23]. Although still nontrivial, the digitalization also allows alignment of serial sections into a common coordinate space, resulting in histology data represented in 3D according to the original shape and structure of the sample. Workflows for the alignment, also known as image registration, and reconstruction from the serially sectioned samples, as well as for visualization of the tissue in 3D are now needed to understand the histological 3D information [11], [14], [17], [24], [25].

Histopathological evaluation of samples is traditionally a very subjective method, where decisions are based on visual inspection of the tissue. To support visual assessment, tools enabling objective recognition and quantification of pathological changes in tissue histology are desired. Computing and visualizing quantitative histological features would also increase our ability to combine histology with other types of data, e.g. information on tumor-associated mutations or gene expression obtained by sequencing, and to explore the effects of tumor genotypes and molecular alterations with respect to other properties within tissue. To be able to quantify changes in tissue morphology, measurable descriptors of the pathology in question need to be computed from imaging data. Numerous approaches for feature-based analysis of histological images have been presented in the literature, see [2], [5], [6], [9], [13], [34], [44] for select examples of earlier studies. In our prior work, we have shown how feature engineering and machine learning enables prediction of neoplastic lesion aggressiveness [29] and in [40] we used a large generic feature set for analysis of spatial heterogeneity and for predicting both the genetic background and spatial location of early tumors in 2D. Furthermore, we expanded this feature set to handle metastatic tissue detection from human breast cancer lymph node samples [38] and further expanded the framework to jointly use both engineered and deep neural network derived features [39]. However, these methods are dependent on the sampling location of the sections and typically do not consider the 3D spatial context in any way.

Native volumetric imaging modalities, such as X-ray microtomography [41], enable straightforward inclusion of 3D spatial context, but lack the resolution for morphological separation of standard histology accompanied by high-resolution imaging. Also with histology, different microscopy modalities can be used for improved specificity, such as by using immunohistochemical staining for labeling specific structures in tissue. With multiplexed staining, advanced visualization and spatial analysis of cellular organization and tissue microenvironment can be done using standard fluorescence image analytics [31], [35]. Modern non-invasive techniques, such as tissue clearing [10], [26], enable harnessing the potential of 3D pathology [18] and enable developing advanced analysis pipelines in 3D [3], [8], but unlike standard sectioning-based histology, such techniques are not routinely available.

In routine histology pipeline, however, specific immunostainings or fluorescence stainings are not always available. On the other hand, using hematoylin and eosin (HE) based standard reference staining of sections only, both the task of capturing the 3D information through registration of sections and the subsequent image-based quantification of spatial organization and tissue morphology become significantly more challenging. Studies by Ourselin et al. [22], Wenzel et al. [42], Arganda et al. [1], Roberts et al. [27], Rosetti et al. [28], Kiemen et al. [12], and [43] show the potential of computational approach in histology-based analysis of the 3D morphology and spatial tissue environment. Methods for detailed quantitative characterization and spatial analysis of HE-stained tissue in 3D are still needed.

Contribution of the current study

The aim here is to develop a method for computational reconstruction, visualization and detailed quantitative assessment of histopathological alterations in three dimensions. We build on the 3D reconstruction of serial sectioned tissue samples presented in our earlier work [11], where several registration algorithms were evaluated and compared. Here, we develop an improved, generic registration framework using the Elastic Stack Alignment algorithm [30] and Bayesian optimization based warm start for hyperparameter tuning from [11] as a basis. Further, we extend the spatial characterization of tissue in 2D which we developed for predicting pathologies and genetic background from murine lesions [29], [40] as well as for human breast cancer metastasis detection [38], by appending the feature set with 3D morphology descriptors and by utilizing full 3D tissue data instead of individual 2D sections for quantification of the intensity and texture features.

Here, we set out to combine inspecting the gross anatomy of tissue and its regions of interest with studying the detailed histology of the samples in quantitative fashion in 3D. We present a series of computed 3D features capable of describing a plethora of tumor characteristics in tissue, as well as present ways to visualize histological features in three-dimensional spatial maps on organ and on a specific region of interest level, providing a novel, detailed view on the spatial properties of 3D tissue samples. We present a modular computational pipeline covering all analysis phases from serial sectioned digital pathology samples digitized as WSIs into organ-level 3D reconstruction of the tissue histology, and quantitative analysis and visualization of tissue morphology and characteristics while preserving the spatial context. We provide source codes and 3D reconstructions for use and to facilitate further development by the community. Using prostate cancer models as a case study, we show that the methods presented here allow inspection and identification of quantitative characteristics relevant to 3D properties of tumors, and identify spatial nuclear gradients in prostate cancer. Although presented here for a specific cancer type, the methods are fully generalizable to different tissue types and pathologies.

Materials

Experimental model and subject details

The data represent prostate tissue samples from FVB/N mice heterozygous for tumor suppressor Pten and either expressing or not expressing transgenic miR-32 (Pten [4] and Pten+/-xARR2PB-miR32 [15]). All animal experimentation and care procedures were carried out in accordance with guidelines and regulations of the national Animal Experiment Board of Finland, and were approved by the board of laboratory animal work of the State Provincial Offices of South Finland (licence number ESAVI/6271/04.10.03/2011). The samples have been previously described in [15], [40].

Methods

Tissue images and ROI separation

The imaging with whole slide scanning device has been previously described in [15], [40]. In brief, the prostates were fixed in PAXgene™ tissue fixative and embedded in paraffin. 5 μm tissue sections were cut, attached to glass slides, and every tenth section was HE-stained. HE-stained slides were whole slide imaged with Zeiss Axioskop40 microscope (Carl Zeiss MicroImaging, NY, USA) with 10x, 20x or 40x objective and a CCD color camera (QICAM Fast; QImaging, Canada) and a motorized specimen stage (Märzhäuser Wetzlar GmbH, Germany). The automated image acquisition was controlled by the Surveyor imaging system (Objective Imaging, UK). Uncompressed bitmap output was converted by JVSdicom Compressor application to JPEG2000 WSI format [36].

Regions of interest (ROI) were manually marked using a freehand selection tool in ImageJ software (National Institutes of Health, Bethesda, MD, USA) [33]. The resulting binary mask was used for extracting the ROI from the full resolution original HE image for further processing. Where one tumor could reach several acinar lumen within a certain section, all affected lumen were included in a single mask. Snapshot images for 2D histology figures were obtained through JVSView program [36], [37] and ImageJ.

Image registration and 3D volumes

Each image stack has approximately 30 sections with 50 μm distance. The remaining tissue sections were reserved for subsequent immunohistochemistry stainings and for molecular measurements. The data was pre-processed essentially as in Kartasalo et al. [11]. In brief, all images were brought to scale of 20x magnification. Sections were segmented from each slide into separate image files by computing Laplacian of subsampled images (1/16 of original dimensions) and thresholding the result using Otsu's method [21]. Spurious objects and holes smaller than were removed using morphological operations. The images were subsampled to /pixel (7849x6529 ± 2055x1759 pixels) corresponding roughly to 10x magnification and saved in lossless tiff format due to limitations in the registration method implementation. The implementation can only produce outputs smaller than pixels. Despite this limitation, the resolution is sufficient for detailed quantitative charaterization of the tissue. Following the preprocessing, images were registered using Elastic Stack Alignment (ESA) [30] using the Fiji [32]/ImageJ implementation. We opted to use ESA in this study based on theresults of a rigorous quantitative evaluation conducted in out earlier work in [11]. Specifically, ESA was found to be the only method among the compared algorithms to consistently provide top scores across majority of the metrics. ESA is an elastic i.e. deformable registration method that applies piecewise affine transformations allowing the formation of complex locally varying transformations. The algorithm first makes a rough initial alignment by computing SIFT features [19] and then finds rigid transformation parameters with RANSAC [7]. Once histological sections are roughly registered, corresponding locations between sections are searched more accurately using normalized cross-correlation. These matches are processed through multiple filters to reduce the number of ambiguous matches. Remaining matches are used to construct a simulated spring system of the registered images, where matching locations across sections are joined together with virtual springs. Such springs are also attached within each section in a triangular mesh pattern to model the tissue's physical properties. The final nonlinear registration is achieved by relaxing this spring system so that the sum of spring forces at every spring intersection is zero. The final transformation is obtained from the triangular spring meshes, where every triangle denotes a piecewise affine transform. (See Supplementary figure 1).

In total, the ESA algorithm has over 30 hyperparameters, posing a complex optimization challenge. Fortunately, its parameters have the tendency to exist in plateaus where the majority of parameter values yield good results when used with histological data with similar characteristics as in our current study [11]. We used a Bayesian optimization based warm-start for parameter optimization, initializing the parameter settings based on the optimum from our earlier study where similar samples were used [11], and then evaluated a fixed amount (N=71) of parameter values around this plateau. We evaluated the result for each combination to select the best registration according to our evaluation metric where

Here, and are consecutive sections in the registered image stack and L is the Laplacian operator. The approach emphasizes tissue area and its prominent features giving minimal weight to the background while not requiring any landmarks. The applicability of the metric was confirmed by a histology expert who visually inspected registered image stacks together with their evaluation scores.

Tissue was segmented from the background to produce tissue masks. Since both the ROI masks and tissue masks were created for the unregistered images, the same registration transforms were applied to the mask images to bring them into the common coordinate system. The tissue masks are needed for efficiency to compute histological features only from tissue regions and also for computing distance features from prostate borders.

Computation of 3D features

From the acquired 3D histological volumes, multiple different types of features were computed. We used a set of histological features modified from earlier studies [29], [38], [40] as a starting point, and tailored the feature set for 3D analysis and extended this set with shape and distance features extracted from 3D volumes of ROI. The features are presented in Table 1 with brief descriptions. The purpose of the 3D feature set is to provide a multidimensional quantitative characterization of the volumetric appearance of the tissue (here: prostate) both in terms of shape and in terms of spatial distribution of ROI (tumors) with respect to the organ and other tumors.

Table 1

Feature descriptions.

Feature	Description
volume	voxel volume * number of white pixels in masked tumor region ( μm^3)
surface_area	average voxel face area * number of voxels in edge volume ( μm^2) (Edge volume is computed by subtracting binary eroded tumor region mask from the original mask and stacking resulting images into a volume.)
sphericity	volume / surface area
dist_section_center_adj_sum	the distance traveled when traversing through each section's masked tumor region's center of mass
dist_section_center_endpoints	the distance between first and last sections' masked tumor regions
dist_section_center_straightness	the distance between first and last sections masked tumor region / the distance traveled when traversing through each section's masked tumor region's center of mass
dist_section_center_*	statistics computed from distances between adjacent sections' masked tumor region centers
dist_section_center_diff_*	statistics computed from differences of distances between adjacent sections' masked tumor region centers
length_pca1	length of the volume along principal axis 1
length_pca2	length of the volume along principal axis 2
length_pca3	length of the volume along principal axis 3
length_pca_21_ratio	length of principal axis 2 / length of principal axis 1
length_pca_31_ratio	length of principal axis 3 / length of principal axis 1
length_pca_32_ratio	length of principal axis 3 / length of principal axis 2
length_pca1_scaled	length of principal axis 1 / length of principal axis 1
length_pca2_scaled	length of principal axis 2 / length of principal axis 1
length_pca3_scaled	length of principal axis 3 / length of principal axis 1
dist_pca1_axis_*	statistics of distances from each voxel to nearest principal axis point 1
dist_pca2_axis_*	statistics of distances from each voxel to nearest principal axis point 2
dist_pca3_axis_*	statistics of distances from each voxel to nearest principal axis point 3
pca1_moment_of_intertia	volume's moment of inertia w.r.t. principal axis 1
pca2_moment_of_intertia	volume's moment of inertia w.r.t. principal axis 2
pca3_moment_of_intertia	volume's moment of inertia w.r.t. principal axis 3
center_of_mass_z_in_prostate	z coordinate of lesion's center of mass in prostate's coordinate system
center_of_mass_y_in_prostate	y coordinate of lesion's center of mass in prostate's coordinate system
center_of_mass_x_in_prostate	x coordinate of lesion's center of mass in prostate's coordinate system
moment_of_inertia	moment of inertia / volume's mass
dist_to_tumor_center_*	statistics computed from distances of voxels to tumor's center of mass
bounding_cube_volume	volume of the bounding cube
bounding_cube_diagonal_lentgh	length of the diagonal of bounding cube
bounding_cube_x_to_diag_ratio	bounding cube x / diagonal length
bounding_cube_y_to_diag_ratio	bounding cube y / diagonal length
bounding_cube_z_to_diag_ratio	bounding cube z / diagonal length
bounding_cube_dim_*	statistics of bounding cube dimensions
convex_hull_surface_area	surface area of convex hull
convex_hull_area_ratio	surface area / convex hull surface area
convex_hull_volume	convex hull volume
solidity	volume / convex hull volume
section_perimeter_sum	sum of sections' masked tumor area perimeters
section_perimeter_*	statistics of sections' masked tumor area perimeters
section_perimeter_adj_diff_*	statistics of differences between adjacent sections' masked tumor area perimeters
dist_nearest_tumor	distance to nearest tumor
dist_furthest_tumor	distance to furthest tumor
dist_tumor_average	average distance to other tumors
dist_tumor_std	standard deviation of distances to other tumors
dist_prostate_anat_center_to_tumor_com	distance from prostate's anatomical center to tumor's center of mass
dist_prostate_anat_center_to_tumor_border	distance from prostate's anatomical center to tumor's nearest border
dist_prostate_border_to_tumor_border	shortest distance between prostate's border and tumor's border

Visualization of features in 3D

The registered image stacks were visualized in Matlab (9.6.0.1114505 (R2019a) Update 2). Slices were drawn along the z dimension, and the same was repeated for x and y dimensions using virtual slices to produce a volume viewable from any angle. Further, the remaining background in the sections was excluded from visualization by setting the background completely transparent while the rest of the numeric range was mapped linearly to make brighter voxels in the volume more transparent. Volumetric visualizations with various effects (movement/rotation of volumes) were then constructed for easy exploration of 3D data in organ and tumor levels, and both using original tissue images and quantitative features as the visualized variable.

Quantification and statistical analysis

The quantification and numerical analysis was conducted using Python software (Python 3.6.8, numpy 1.16.0, pandas 0.23.4, scikit-image 0.14.2, all details in environment.yml in the Supplementary code repository), and visualizations in Matlab (9.6.0.1114505 (R2019a) Update 2) as described Method Details. Description of the computational experiments are given in Results. The presented data consists of prostates of 6 mice, with 2-6 tumors per prostate. Half of the mice were wild type and half transgenic for overexpressing miR-32 in their prostates. The sample size of 22 tumors spans a representative spatial distribution in prostate.

Results

Creating 3D representations of mouse prostates with tumors from histological sections

We set out to create a pipeline for quantitative spatial analysis in 3D (Fig. 1). Mice heterozygous for tumor suppressor Pten form mouse prostatic intraepithelial neoplasia (mPIN) within 8–12 months (Di Cristofano et al. 1998). In here, we used prostate samples from Pten mice of 10-11 months, when recognizable tumors are evident. Our dataset included images of prostates from 6 mice, harboring between 3-9 HGPIN tumors per prostate [15], [40]. Three of the mice were transgenic for overexpressing miR-32 in their prostates [15]. The distribution of these tumors in the prostate has previously been shown to be mostly located to the lateral prostate (69%), with individual tumors in the dorsal (22%) and ventral (9%) prostates. To ensure sufficient volumes and 3D structure to be analyzed, only tumors spanning at least four consecutive images ( distance) were included in the 3D analysis, resulting in 22 tumors in the dataset (2-6 tumors per prostate).

Figure 1

Pipeline of quantitative spatial analysis in 3D. A) Serial section WSIs are aligned using elastic registration into a 3D stack. B) The 3D tissue stack is used as a whole organ model, and regions of interest (tumors) are extracted for subsequent analysis. Binary masks corresponding to the whole tissue and ROI areas are collected for visualization and analysis purposes. C) Quantitative analysis is performed for the 3D volumes using shape analysis of the ROI volumes and spatial histology analysis for the 2D images within the stack. The quantitative features can be explored using multidimensional numerical analysis and by visualizations in their 3D spatial context.

To visualize the histological information of the tissues in three dimensions, we first performed registration of WSIs to image stacks by applying an automated elastic registration pipeline using ESA as the registration algorithm, with parameter settings optimized for the task. The six prostates used in this study are visualized as histological image stacks in Fig. 2 (left panels), showing the appearance of the whole tissue as 3D objects and visualizing the main characteristics of the organ. To better detect the objects of interest in the organs, transparent visualizations of the organs are built, with tumors indicated with individual colors (Fig. 2, right panels). These visualizations provide an easy and intuitive way to study the overall context of the tissue as well as individual tumors for their locations, shapes and sizes. Especially, rotating the volumes to enhance the visibility through observation via multiple angles enhances examining the tissue and tumor characteristics (Supplemental video 1, available on supplementary site).

Figure 2

3D reconstructions of prostate tissues. Six mouse prostates used in the study are shown. The prostates represent Pten+/− (prostates 1-3) or Pten+/-xARR2PB-miR32 (prostates 4-6) genotypes in FVB/N mouse strain. Left: Histological 3D reconstructions composed of HE-stained sections ( at intervals) For visualization, standard surface interpolation is used in spaces between sections. The prostates are presented the urethra (white asterisks) with its prominent muscle layer facing up, with ventral prostates in the front right corner. Due to the unencapsulated nature of mouse prostate, the gland structures are readily visible whenever not masked with adipose tissue (black asterisks). Right: 3D constructions showing the outer borders of the tissues in sections, as well as the tumors with individual coloring based on 2D tumor masked regions of interest for each prostate.

Quantitative 3D features in tumors

We have previously performed quantitative analysis of the histology of these tumors from the annotated tumor areas in 2D section WSI's [40]. While this analysis provided valuable information of the histological characteristics of the tissues and tumors, information of the 3D context, shape, and spatial and volumetric characteristics of the tumors could not be considered with such an approach. Here, we wanted to quantitatively characterize the three-dimensional aspects of these tumors. For this, we engineered a set of features to describe the 3D nature in terms of sizes, shapes and different aspects of variance in search for quantitative determinants describing and possibly capable of subgrouping the tumors (Supplementary Table 1). In addition, we calculated statistical descriptors for certain features, such as mean, median, minimum, maximum, sum, variance, standard deviation, kurtosis, and skewness, to capture the distribution of the feature values.

Further, in addition to basic size-describing features, such as volume and features derived from dimensions of the tumors (such as tumor lengths, widths, and thicknesses calculated from projections on principal axes), we included measures describing and helping to identify types and variance in tumor shapes. These features are derived from surface area, bounding cube and convex hull, such as sphericity and solidity. Also features utilizing aspects of the area and number of sections occupied by the tumors, such as features derived from the tumor sections' (i.e. tumor areas in 2D in each tissue section) perimeters were included. To assess the locations of the tumors in relation to the host tissue as well as to each other, we engineered distance features using the anatomical center point in each tissue as a reference point. These anatomical center points, defined manually by an expert, are visualized in Supplemental Figure 2, depicting also the tumor locations in the prostates with views from lateral and dorsal sides.

3D features distinguish tumor subgroups

We computed the values for the 3D features for the 22 tumors in the dataset and examined their capability to distinguish tumors into subgroups. Clustering analysis (Fig. 3A) shows that the main feature cluster is formed by approximately half of the features more or less correlating with the volume of the tumors, while the second half of the features forms smaller and less dense clusters. Clustering based on tumors, on the other hand, forms three major groups directed by several features describing the size and shape of the tumors (Fig. 3A). We used principal component analysis (PCA) to reduce dimensionality of the data in order to determine whether the features reveal subgroups relevant for the biology of the tumors. A PCA plot of individual tumors based on the 3D features shows that the tumors in the cluster groups indicated in Fig. 3A separate to distinguishable locations (Fig. 3B). Closer inspection of 20 features with the highest weight in the PCA analysis (Fig. 3C) indicates that cluster 1 exhibits, for example, relatively high and varied bounding cube dimensions, high voxel distances to the first principal axis, high length of the third principal axis, high distances and variance in tumor area centers between sections and low solidity. These attributes indicate large but also complex shaped tumors. On the other hand, tumors in cluster 2 have relatively low distances between tumor section centers, long principal axes 2 and 3 as well as high ratios between the lengths of these axes and principal axis 1, indicating more round tumor shapes. The rest of the tumors, corresponding to cluster 3, are relatively more compact, exhibiting smaller size according to several features and high solidity.

Figure 3

Computational 3D features group prostate tumors. 3D features were computed for each tumor. A) A heatmap representing 3D feature values for individual tumors analyzed. Clustering of the data reveals three distinctive groups of tumors, marked as clusters 1-3. B) Principal component analysis (PCA) plot of 3D features indicating individual tumors. Clusters 1 and 2 indicated in A are distinguished. C) Biplot of the PCA plot in B, showing 20 features with most weight in positioning the tumors in the PCA. D) The same PCA analysis as in B, indicating the genotypes of the tumors. Interestingly, while the clustering in A does not indicate genotype-dependence, large tumors in the dataset are positioned to distinguishable groups according to their genotype as indicated.

Although these Pten tumors belong to two genetic groups (miR-32 wt and miR-32 transgenic/overexpressing) and show certain prostate cancer marker expression differences, they are indistinguishable by visual inspection in conventional 2D histopathology [15]. Here, inspecting how tumors of the two different genotypes are positioned in the PCA analysis revealed that, while small tumors of both phenotypes are intermixed, the 3D feature values show separation among the genotypes in the large tumors (Fig. 3D). The larger tumors likely represent an advanced growth phase compared to the small tumors corresponding to early tumors or tumors with less growth potential. Thus, these results suggest that the tumor growth patterns have distinguishable features in 3D analysis that can be linked to their genotype, depending on the growth phase.

The differences indicated in the clustering and PCA analysis prompted us to investigate more thoroughly selected 3D features in relation to the tumor phenotypes (Fig. 4). Especially, the size and shape-related features represented strongly in the clustergram heatmap analysis, as well as certain shape- and location-related features estimated to be useful spatial descriptors for tumor pathology, were investigated more closely. Three basic 3D object features (volume, sphericity, and solidity) are useful in providing individual descriptive patterns for each tumor. Interestingly, sphericity is increased overall with increasing volume of tumors (Fig. 4B), indicating that larger tumors tend to become more round rather than having lengthy shapes. On the other hand, tumor solidity is decreased with increasing volume, although in a nonlinear fashion (Fig. 4C). These notions support the idea that small and early in situ tumors grow along the free space in the narrow and winding prostate gland lumen, whereas the increased growth pressure makes way in the tissue as the tumors enlarge and get more complex surface shapes. To this end and as already indicated in the PCA in Fig. 3B, using the principal axes of tumor maximal lengths in the first, second and third angle, in addition to their correlations, create useful features describing and distinguishing the shapes of the tumors (Supplementary Figure 3).

Figure 4

Tumor level inspection of selected size, shape and location features. A) Selected 3D features representing size, shape and location of tumors in the prostate features, and volumetric surface visualizations of each tumor from two different angles. B)-G) Correlation plots between selected size, shape and location features. See main text for details.

The features describing the tumor locations in the organ show heterogeneity with respect to inter-tumor distances (Fig. 4A,D). Also tumor distances to the anatomical center and prostate edge borders vary (Fig. 4A), as expected in an organ containing separate lobes with variable length. This may explain why the distance to the anatomical center is not correlating with the above mentioned features (Fig. 4A,E), while the distances to edges of the organ are clearly shorter with large tumors (Fig. 4F,G). This may be indicative of the fact that the prostate in mouse is an unencapsulated organ, in which increasing spatial pressure due to a growing tumor is most easily relieved by directing growth towards the abdominal cavity.

Visualization of quantitative histological analysis in 3D

Many histological entities, such as pathological areas and tumors, span more than one histological 2D section, and the histology may vary between different areas of the region of interest. We have previously performed a computational analysis of 2D histological features of the tissue sections used in this study. Here, in order to understand histological features in their 3D tissue context, the 3D representations of the tissues with tumors are supplemented with this quantitative 2D histological information.

We selected two differing, informative 2D features from our previous study [40] to visualize in the 3D representations and study the distribution of quantitative data in the whole prostates and the tumors. Eosin channel intensity quantifies the red signal (eosin signal computationally separated from hematoxylin and mean signal value quantified within local support area) highly enriched in proteins in the cell cytoplasm, collagen, and muscle fibers. Nuclear density is based on nucleus detection from hematoxylin signal computationally separated from eosin and performed within local support area, and density calculated through windowed kernel density estimation.

Fig. 5A shows these features as 3D representations of the same prostate, indicating the differential distribution of the signals along the tissue and visualizing different structures. While the eosin channel intensity distinguishes the urethral muscle wall especially well (Fig. 5A,B,level A), it is enhanced also in tumor tissue (Fig. 5A,B, level B). Nuclear density highlights especially the intraurethral glands from normal tissue (Fig. 5A,B,level A), and this feature is also well enhanced in the tumor tissue (Fig. 5A,B, level B). Both features visualize the tumor already in the 3D representations, more clearly by the nuclear density feature (all 6 prostates shown in Supplementary Figure 4).

Figure 5

Quantitative histological features computed from 2D histology visualized in 3D format. A) 3D representations of a prostate from 1) HE-stained sections showing tumors in color and two example levels of histological sections (Level A and B, upper panel), 2) eosin channel intensity feature (lower left panel), and 3) nuclear density feature (lower right panel). B) Examples of histological HE-stained sections (left panels) and intensity plots of eosin channel intensity and nuclear density features (right panels) of the sections indicated in A). Upper two rows show sections of a whole organ, lower two rows show magnification of inserts marked in the upper sections. In level A, urethral muscle wall with high intensity in eosin channel (black asterisk) and intraurethral glands with high intensity in nuclear density feature (while asterisk) are shown. In level B, a tumor area with distinctive density compared to normal glandular environment in both feature channels is shown (arrowhead).

A closer inspection of the feature representation in 3D focusing only on the ROI may allow observations not visible with surrounding tissue masking the view. Studying the nuclear density feature in tumors only (Fig. 6 and Supplementary Figure 5) reveals that nuclear distribution is not even throughout the tumors, but that there are distinct density gradients. Observing the tumors with the feature from several angles and comparing that to 3D representation of the tissue with HE-staining reveals that the densities of nuclei are higher closer to the outer edges of the organ. This notion indicates that the growth density of the tumors is higher towards the edges of the unencapsulated organ and, thus, closer to lower pressure by the surrounding tissue.

Figure 6

Prostate tumors have nuclear density gradients revealed by 3D visualization. An example prostate (no 6) visualized in 3D based on A) computationally determined nuclear density and B) the corresponding HE-stained sections with indicated tumors. C-E) Tumor areas are visualized for their nuclear density. The magnified views and angles in D-E) reveal a gradient-like spatial density alteration of nuclei in a large tumor, with denser areas with higher number of nuclei close to organ edge and likely corresponding to areas of increased growth rate.

Discussion and conclusions

Computational methods for visualization and quantitative assessment of histopathological alterations in three dimensions are required to better understand pathological changes in the tissue environment. In here, we reconstructed 3D representations of whole organs from serial histological sections and analyzed the tissue characteristics and regions of interest in 3D. Our work shows that quantitative assessment of tumor sizes, shapes, and separation between different spatial locations within the organ can be used to assess tumor characteristics, and to group tumors according to their 3D properties. In addition, our presented approach for reconstructing the three-dimensional structure of tissue, extracting quantitative descriptors capturing its spatio-morphological characteristics, and for visualizing both the histology and computationally quantified feature representation in 3D, provides an intuitive way to observe pathological characteristics. Particularly, it offers enhanced view on an organ or, e.g., tumor level exploration of the studied tissue, as opposed to the traditional way of inspecting the sample section by section in 2D.

Traditionally, quantitative features have been used in machine learning algorithms, and many methods enable defining which features have the highest influence on classification results. We have previously used this prostate tumor material in a study where we analyzed spatial heterogeneity of neoplastic alterations in these prostate tumors in a 2D setting [40]. We computationally quantified hundreds of features from the histological WSIs and used them to build a feature-based machine learning model to separate genetic alterations between tumors. Now, we presented additional quantitative features to describe the tissue and its alterations in 3D. Furthermore, we presented 3D visualizations to explore the spatial distribution of the important features to associate them to anatomical locations and to study what properties they capture. The 3D visualization and exploration method could be supplemented with additional types of data, for example, immunohistochemical stainings, spatial transcriptomics, and volumetric imaging modalities.

We used mouse prostates with prostate cancer tumors as a use-case. With 3D representations, the tumors could be explored in ways that are not possible through standard 2D visualization of histological sections – not even when several sections are utilized in a row. The ability to observe the tissue and the tumors from multiple angles clearly enhances the information gained of tumor distribution to tissue, their growth patterns and their heterogeneity. Our work presents methods and descriptors to study and represent tumor shapes and sizes, information which, in here, can be used to distinguish between genetically different groups of tumors. Combined to increased number of samples and large-scale molecular information, our methods for quantitative characterization of tissue in its 3D spatial context harbor vast potential to investigate and understand the effects of e.g. different genetic alterations to tumor biology.

Our results underline that there is heterogeneity in composition and cellular organization within individual tumors, even in the relatively homogeneous genetic setting utilized in here. Thus, the quantitative computational methods, such as those presented here, are required to screen for relevant descriptors out of reach to the human eye. In our example case, the two different genetic groups of mice (Pten and Pten+/-xARR2PB-miR32) showed indistinguishable histopathology [15], while the computational 3D feature analysis here revealed differences between the larger tumors in the two genetic groups. Furthermore, our visualization of 2D histological features in 3D showed increased utility over quantitation only, revealing that prostate tumors have nuclear density gradients indicating areas of tumor growth directions. These valuable biological insights provide novel information to tumor biology in tissue, and serve as points of direction for future studies, especially considering that the tumors in this Pten prostate cancer model are multifocal, representing well the situation in human prostate cancer.

Our work provides proof-of-principle for gaining significant additional information from histology by studying it quantitatively in 3D. We used mouse prostates in this study, but the methods presented here are applicable also to human tissues as well as several other types of pathologies in addition to cancer. We envision the primary use of 3D reconstruction from histology as a very promising and powerful technique for answering questions on spatial pathological patterns in tissue and for characterizing spatial heterogeneity of the tissue environment. Detailed and quantitative spatial characterization of tissue environment also enables the use of advanced visualization techniques, such as virtual reality based immersive and interactive exploration of feature data in 3D [16].

Incorporating various other measurement modalities in the 3D models alongside the histology-based examination shown here, can provide detailed spatial models combining histology with, e.g., molecular level measurements. For example, multiple immunostainings and spatial transcriptomics data would provide detailed spatial models with molecular level characteristics of the studied samples. Research applications aiming at detailed spatial mapping of tissue environment, such as characterization of tumor growth patterns, form the most prominent application areas, whereas clinical use, despite the increased availability of enhanced automated pipelines for histological processing and routine slide scanning, is limited by the requirement for serial sectioned sample volumes. However, in contrast to methods requiring intact tissue for volumetric imaging, 3D reconstruction and analysis based on serial histological sections can also be utilized for samples already processed for histology, or e.g. archived clinical samples. This opens up enormous possibilities and sample numbers to utilize in further research.

Limitations of the study

Despite presenting a generic pipeline for quantitative 3D characterization of histological samples, here we limited our analysis to mouse prostate tissue in an organ-level case study. Applicability in studying pathologies in other tissues and diseases requires 1) sufficient number of serial sections to cover histologically meaningful volume enabling extraction of quantitative features in 3D, and 2) sufficient computational capacity to handle the memory and computationally intensive pairwise registration and computational feature extraction phases. The implementation (see: Data and Code Availability) is modular and requires a versatile computational software stack. Our study covers exploratory analysis and interpretation of the feature data and volumetric visualizations, while, e.g., statistical significance based feature selection or correlative analysis of genotype-feature (phenotype) associations need to be explored in subsequent studies.

Declarations

Author contribution statement

Pekka Ruusuvuori, Leena Latonen: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Masi Valkonen: Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Kimmo Kartasalo, Mira Valkonen: Contributed reagents, materials, analysis tools or data.

Tapio Visakorpi, Matti Nykter: Conceived and designed the experiments; Contributed reagents, materials, analysis tools or data.

Funding statement

This work was supported by the , projects #313921 & #326463 and #314558 & #326364 (PR), #279270 & #317755 (TV), #312043 (MN) and #317871 & #334774 (LL), (PR, TV, LL), ERA PerMed 2019-2022 - ABCAP project (PR, LL), (TV, LL), and Foundation of the (LL), and (UTUGS) and Turku University Foundation.

Data availability statement

Data associated with this study has been deposited at https://gitlab.com/BioimageInformaticsGroup/3d-histology/.

Materials availability

This study did not generate new unique reagents, mouse lines or tissue samples.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.