Elastic transformation of histological slices allows precise co-registration with microCT data sets for a refined virtual histology approach

Although X-ray based 3D virtual histology is an emerging tool for the analysis of biological tissue, it falls short in terms of specificity when compared to conventional histology. Thus, the aim was to establish a novel approach that combines 3D information provided by microCT with high specificity that only (immuno-)histochemistry can offer. For this purpose, we developed a software frontend, which utilises an elastic transformation technique to accurately co-register various histological and immunohistochemical stainings with free propagation phase contrast synchrotron radiation microCT. We demonstrate that the precision of the overlay of both imaging modalities is significantly improved by performing our elastic registration workflow, as evidenced by calculation of the displacement index. To illustrate the need for an elastic co-registration approach we examined specimens from a mouse model of breast cancer with injected metal-based nanoparticles. Using the elastic transformation pipeline, we were able to co-localise the nanoparticles to specifically stained cells or tissue structures into their three-dimensional anatomical context. Additionally, we performed a semi-automated tissue structure and cell classification. This workflow provides new insights on histopathological analysis by combining CT specific three-dimensional information with cell/tissue specific information provided by classical histology.

Conventional histology of paraffin embedded tissue specimens has been the gold standard for the analysis of tissue samples for decades. While countless staining protocols and the use of microscopes allow very specific analysis down to (sub-)cellular level, histology falls short when it comes to acquiring 3D information of the tissue. Serial sectioning marginally circumvents this problem, but it is very time consuming, labour intensive and fails to deliver a real 3D dataset 1,2 . Therefore, techniques that can correlate cellular and morphological information within the 3D context of the entire tissue volume are in great demand.
MicroCT imaging can deliver high resolution 3D datasets 3 . Soft tissue specimens, however, require additional contrast commonly achieved using heavy ion-based staining methods. Our group, among others, has used different CT staining protocols to improve the specificity of microCT imaging. While phosphotungstic acid or iodine-based protocols can greatly increase the CT contrast of tissue, and therefore allow for a 3D assessment on a micrometer scale 4-8 , they do not increase the specificity to identify certain cells or tissue structures. Other approaches involve using reagents known from classical histological stainings. Eosin can be used to stain cytoplasm 9 , while hematein (led-III-acetate) functions as a nucleus specific staining agent 10,11 . While these approaches seem to take the specificity of microCT based virtual histology to a new level, they are hindered by limited sample sizes as well as by the fact that only a few applicable protocols exist to this day. Metscher et al. showed that a specific antibody staining of chick embryos for microCT analysis is possible using silver labelled antibodies 12 . Due to the limited tissue penetration depth of the metal labelled antibodies OPEN 1

Results
To develop the workflow of the co-registration between label free 3D microCT based virtual histology, classical histology and immunohistochemistry, we analysed tumour tissues from a mouse model of breast cancer by SRµCT and subsequent histological processing.
Workflow for precise registration of histological sections and SRµCT data with Fuxlastix. Initial overlay of the CT and histological images required finding the position and orientation in the CT data set, which corresponds to the actual paraffin section, cut by the microtome. To achieve this, the CT data set was virtually cut with a plane that was manually adjusted until the virtual slice matched the histological slice. In order to compensate for the deformation due to cutting and staining procedures, a tailored pipeline was developed, as illustrated by the flowchart shown in Fig. 1. The virtual histology slice obtained from the 3D CT data set was used as a 'fixed image' template, representing the correct morphology before the histological processing. The histology slice was considered the 'moving image' , representing non-uniform deformation upon the completion of histological processing. We then used a b-spline interpolation, which represents a grid that allows stretching and moving each position individually. Although useful for non-homogenous deformation, this model cannot handle strong global mismatches, especially rotations, as they are always approximated by deforming the grid. Therefore, initially the translation and rotation of the moving image with respect to the fixed image were corrected using a Fourier-Mellin algorithm 30,31 . For simplicity, either a single colour channel or the average of all channels from the moving image was used further. In addition, the histology images were inverted so the relevant image content is given as high grey values for both histology and CT. Both, the moving and the fixed images, were downscaled to 1/10th of their original size before registration with Elastix was performed. Mutual information 32 was used as optimisation criteria and a b-spline as the deformation model 33 . The mutual information used as metric for registration is based on pairing the intensities at a certain pixel in both the fixed and moving image. To suppress the influence of noise, the intensity range was binned to 32 bins respectively. The transformed channels were then merged into a single colour image. The resulting deformed image is shown in Fig. 1 (two bottom images). The blue border (Fig. 1, bottom left) indicates the degree of deformation calcu-  www.nature.com/scientificreports/ lated by the registration pipeline. The checkerboard view demonstrates the precision of the registration (Fig. 1, bottom right).
As an example, to demonstrate the capabilities of the pipeline, Fig. 2 illustrates a 3D rendering of a mouse breast tumour that was intratumourally injected with radiopaque barium nanoparticles (BaNPs). The data set was virtually cut open at the position of the subsequently performed histological section. The HE-stained section was registered to the virtual slice and placed as a plane into the 3D architecture of the tumour. The image demonstrates schematically that by using our workflow, specific information, like the occurrence of tumour cells or collagen fibres from the histology can be related to the 3D location within the entire specimen. The BaNPs, which are invisible in histology, can be easily discerned in the CT data set. In addition, 3D features like the position and extent of the "blobs", which represent hollow cavities in the tumour and cannot be effectively analysed by histology, can be easily assessed volumetrically by CT.

Evaluation of the registration quality.
To evaluate if the deformation of the histology resulted in an optimal fit with the virtual cut through the CT data set, we calculated the previously introduced displacement index (DI) 8 . DI assesses the displacement and the mutual information using a block matching approach. An www.nature.com/scientificreports/ ideal match should result in a DI of 0. However, since histology and CT show different image content even for perfectly overlaid images, the DI will never reach this value. Figure 3A shows an overlay of CT data and the green colour channel of the corresponding MTS-stained histological slice (pseudo-coloured in green to improve visibility) before the registration. Evidently, the MTS data was deformed during the cutting and staining process and therefore did not match the CT slice. Figure 3B shows the same overlay after the histology data was transformed using Fuxlastix. It is clearly visible that a near perfect match was achieved. Figure 3C,D illustrate the calculation of the DI of the images pre-and post-transformation with a block size of 500 pixels with an overlap of 50%. The congruence between the CT tiles and the histology tiles is directly proportional to the size of the green circle inside that tile. The red lines show the calculated displacement of each tile. The results for all performed transformations demonstrate a larger DI of 11.7 ± 2.0 (n = 15) in the untransformed data as compared to a DI of 6.9 ± 2.0 (n = 15) in the transformed data sets, confirming a significant improvement of the overlay after registration of the histology data.
Image registration of differently stained serial histological sections with microCT. In order to provide a multi-parametric analysis of tissue samples, different staining protocols are commonly applied on consecutive histological sections. Since every histological slice may display different degrees of non-uniform deformation, a direct overlay provides limited information. To correct for various degrees of deformation between the different histological slices, we used an unaltered CT data set as a fixed image template. We applied our elastic image registration technique on every corresponding virtual cut through the CT data set of mammary tumours (Fig. 4A) and performed separate registration with consecutive slices stained with different techniques. Masson Trichrome staining (MTS, Fig. 4B) was applied to depict collagen fibres in blue. Figure 4E shows a magnification and explicitly demonstrates the continuity of the structures achieved by the registration. Specific structures like the blood vessel were identified through histology and confirmed by tracking the vessel throughout the 3D CT data set. Figure 4C,D show specific IHC stainings for the macrophage and monocyte marker CD68 and for the SV40 T-Antigen (T-Ag), used as specific tumour marker in the WAP-T tumour mouse model, matched to Allocation of NP accumulation sites to anatomical structures in mouse breast tumour samples. In order to evaluate our workflow in the scope of a biomedical question, we applied the developed registration workflow in a NP based theranostic strategy to assess NP distribution as one of the key issues. BaSO 4 containing NPs were injected directly into the mouse breast tumours, followed by low energy irradiation over multiple therapy sessions. The tumours were then explanted, and high-resolution CT and histology were performed as described. The high image quality of the embedded breast tumour in CT was achieved by applying propagation-based imaging without the use of additional contrast agents, which would have reduced the relative contrast of the nanoparticles. BaNPs were easily discerned (Fig. 5A), due to their high X-ray contrast. The BaNPs, however, could not be visualised in histology as shown in the corresponding MTS-stained slice (Fig. 5B). Thus, we applied the same workflow described above in order to place the BaNPs in the anatomical context of the tumour. Figure 5B shows the histological slice corrected for deformations during the cutting process as indicated by the inhomogeneous boundary regions coloured in blue. Figure 5C shows the same histological slice with the overlaid and segmented BaNPs. The BaNPs were segmented by using a simple thresholding method (Otsu's method) and are shown by a yellow pseudo-colour enabling to depict the position of NPs within the tumour.
In the tumour shown here, BaNPs are concentrated in the periphery of necrotic regions, which are characterised by the lack of red staining in MTS. As indicated in the close-up image (Fig. 5D), our registration approach allows to analyse the local environment of small clusters of particles. Since the NPs were designed for contrast enhanced radiation therapy, we assume that the visible necrotic regions were most likely formed due to the increased local radiation damage caused by the particles.
The co-registration shown in Fig. 5E depicts another example of a breast tumour where the BaNPs can be allocated only to the tumour capsule, indicating a different pattern of NP distribution after intratumoural injection. Figure 5F illustrates the corresponding anti-CD68 stained IHC slice after the registration was performed. Furthermore, the fusion with an anti-CD68 IHC stained slice (Fig. 5G) demonstrates that the localisation of the NPs coincides with regions of CD68-positive cells, indicative of macrophages that have phagocyted the BaNPs. Figure 5H shows a magnified detail view of Fig. 5G clearly depicting a co-localisation of the BaNPs and phagocytic cells.
These examples neatly demonstrate how the information unachievable by classical histological analysis can be supplemented by means of registration with high resolution CT.
Precise co-registration of microCT and histology improves automatic cell classification. As an example of how combined analysis of microCT and histology can be used to improve automatic cell or tis- www.nature.com/scientificreports/ sue classification, we used a pixel-wise mixed density model approach. Figure 6 shows a microCT image of a local area of a mammary tumour containing BaNPs (Fig. 6A) and the corresponding part of the MTS-stained histological slice elastically matched using the proposed pipeline (Fig. 6B). Thus, each pixel is characterised by its intensity value of the CT data and its three colour components of the MTS staining generating a 4D feature space. Since it cannot be expected that the registration between the two images is precise on a pixel-wise basis, the feature space showed highly overlapping clusters limiting the automatic classification of different structures. Therefore, we used a 4D Gaussian mixed model with 6 components that was automatically adapted to the data by applying an expectation maximisation algorithm implemented in ski-learn. Consequently, the probability for each pixel to belong to each of the 6 classes is shown in Fig. 6C-H. As a result, the class shown in Fig. 6C represents the tumour cells (25% probability). Extra cellular matrix components, such as collagen are represented by the class shown in Fig. 6E. The fat cells (Fig. 6G) can easily be separated from small vessels and empty areas in the tissue (Fig. 6D), which would not have been possible using only the histological slice (Fig. 6B). Since the fat contained in the adipocytes is washed out during the staining process, the fat cells appear to be indistinguishable from the porous morphology of the tissue, making the ability to segment them even more remarkable. Figure 6F mainly depicts the blue stained fibre content in the MTS staining, while Fig. 6H highlights interfaces and cell membranes between the different cell types. Our results show that a simple tissue classification can be improved by using correlative CT and histology data sets. Without the combined information provided by the CT and the registered MTS histology an automatic classification into the shown classes was not successful, as shown in the supplemental files ( Figure S1, S2).

Discussion
Here we present a novel workflow to precisely register histological sections with virtual cross-sections within microCT data sets of paraffin embedded tissue and demonstrate its performance on tumours of a breast cancer mouse model. We developed a simple to use workflow based on Elastix 34,35 to perform an elastic registration of the data sets, which is available here https:// github. com/ xPITc oding/ Fuxla stix. We show that a near cellular precision between two imaging modalities can be achieved. To simplify the registration, we reduced a 2D-3D problem to a 2D-2D problem by performing a manual virtual slicing of the CT data. While this makes the registration itself much easier, the resulting quality of the overlay is always limited by the precision of the initial virtual sectioning. Even though we believe that an automated co-registration technique would not necessarily offer a more precise match of both modalities, it could speed up the whole workflow significantly. Since cutting of the histological sections is performed by a microtome, which is essentially a static blade, it is reasonable to assume that the histological slice can be described by a 2D plane through the microCT data set and the deformation associated with the cutting is restricted to that plane.
In the literature most methods for image co-registration are either intensity, landmark or segmentation based [36][37][38][39] . Since, in our case, microCT data sets and histology have a vastly different image content a flexible registration method was needed. In addition, both SRµCT and histology are prone to artifacts (air bubbles, paraffin cracks, tissue breakage, uneven staining). Therefore, a very robust registration method was required. Because of this, we used a mutual information-based registration approach, rather than a landmark based method. Chicherova et al. demonstrated successful registration between microCT and histology in 3D using a SURF (Speeded Up Robust Features) algorithm 40 . Since SURF is based on the automatic identification of corresponding landmarks, it appears to work well on histochemical stainings such as the HE staining presented in their publication. In the case of antibody-based staining, in which the image content differs strongly between CT and histology, we speculate that the application of SURF would be challenging. In contrast to SURF's landmarking, our approach uses mutual information and considers the entire image content, which has proven to be more robust and precise in our IHC data. However, Chicherova et al. achieved the best results when a curved plane was fit to the 3D data set, while our approach works under the assumption of a flat plane and 2D registration. Therefore, it is reasonable to further evaluate whether the registration quality can be improved by forfeiting the assumption that the microtome cuts strictly on a plane. Additionally, Chicherova et al. settled on a deformation model based on Legendre-polynomials in order to avoid the flexibility of the b-splines 40 . This is indeed a reasonable approach when applied in dense tissues like bone and cerebellum. However, our specimens displayed a large variability in deformation magnitude at different regions in the tissue. Thus, we believe that the b-spline approach was indeed required in our case to achieve a precise match. On the other hand, we are aware that our elastic registration pipeline cannot be perfectly validated as no ground truth data of an optimal overlay can be generated and elastic image registration will always converge to a solution. In the future, with the availability of more data, the maximal local deformations for a given tissue can be evaluated and used to modify the penalty term to a solution with large local deformations of the b-splines.
It is important to point out, that SRµCT does not interfere with subsequent staining procedures of paraffinized samples commonly employed in clinical routine. SRµCT can deposit very large X-ray doses in the sample (in the range of kGy for the samples shown here). In this study we did not apply any additional staining to increase the X-ray contrast of the sample, therefore the absorption is relatively low, leading to reduced dose deposition. In contrast to that, our former experiments 6 using phosphotungstic acid (PTA) stained tissue samples led to a significantly increased X-ray absorption and therefore dose deposition. SRµCT in PTA stained samples can lead to micro cracks in the paraffin, which can impede the sectioning process. For the samples presented here, no visible crack formation was observed. For both PTA stained and unstained SRµCT scanned specimens no significant impairment of histology and IHC was observed. Paraffin blocks, when compared to other embedding materials, provide a number of advantages such as: (i) access to more staining protocols, (ii) easier handling and (iii) use of subsequent analysis techniques such as atomic force microscopy. Thus, our approach may be integrated into the standard workflow of histopathology.
In BaNP injected mouse tumours we demonstrated that our approach facilitates the correlation of NP location to anatomical features such as fibre structures as components of the extracellular matrix, cell density in necrotic areas and macrophage localisation. However, since histological sections are performed in a serial manner, they do not represent the same cells/tissue structures. The co-localisation to virtual histology slices can therefore not result in an optimal match on cellular level for every slice. New approaches offer the advantage of having multiple antibody stainings on the same slice 41,42 . This would enable an improved multiplexing of microCT data with information from multiple stainings. Therefore, artificial staining of the microCT 3D dataset by using machine learning approaches to transfer the appearance (style) of histology to the microCT data is feasible, even at locations for which a corresponding histological slice does not exist. A similar approach based on style transfer learning has been established for unstained tissue images, where artificial HE and MTS stainings were performed based on the intensity of autofluorescence of the tissue 43 . We are thus confident that a pseudo-histological-stained X-ray based virtual 3D histology could be achieved and would furthermore allow to perform histology evaluation at any position and orientation in no time. Moreover, the specimen will not be destroyed and multiple virtual stainings would be possible.
Since our approach allows a precise match between histology and virtual microCT sections, all information can be assigned to the same structure at near cellular level. Such a location can be described by a multi-dimensional feature vector as demonstrated by the combination of the local grey values in the microCT data with the intensities of the three different colour channels in the MTS-stained slice. These feature vectors enable us to employ unsupervised machine learning techniques to identify different cell and/or tissue types. We demonstrated www.nature.com/scientificreports/ this by the automated separation into healthy tumour regions, fibre regions, fat cells and NP positive regions. The combination of more and different stainings will help to further improve this technique.
In summary, we believe that this workflow based on our open-source frontend for Elastix could aid the histological tissue analysis by supplying additional 3D information, in a large variety of applications in biomedical research. Moreover, it might present a valuable starting point for machine learning strategies and a deeper understanding of pathological processes.

Methods
NanoPET BaSO 4 nanoparticles. Barium sulphate based nanoparticles (BaNPs) were produced by chemical precipitation at ambient conditions 44 . After a purification step, the particles were sterically stabilised using a biocompatible polymer. The obtained highly stable colloidal suspension was sterilised and formulated for in vitro/in vivo use. The hydrodynamic diameter of the particles was around 120 nm.
When the tumours reached a size of 300 mm 3 , 15 µl of BaNPs were injected intratumourally by performing three injections at different sites with a 10 µl glass syringe (Hamilton), followed by irradiation of the tumours with a total dose of 30 Gy over 15 therapy sessions in a small animal in vivo microCT (QuantumFX, Perkin Elmer).
The mice were sacrificed using a CO 2 overdose in combination with a cervical dislocation. The tumours were explanted, fixed in 4% formalin, and embedded in paraffin for further histological analysis.
All animal in-vivo procedures were performed in compliance with the guidelines of the European Directive (2010/63/EU) and the German animal ethics regulations and were approved by the local ethics office (Niedersächsisches Landesamt für Verbraucherschutz und Lebensmittelsicherheit, LAVES, ethics approval 33.9-42502-04-18/3022. Furthermore, the study was carried out in compliance with the ARRIVE guidelines (https:// arriv eguid elines. org/).

Synchrotron radiation microCT (SRµCT).
All high-resolution SRµCT acquisitions were performed using the white/pink beam setup of the SYRMEP (SYnchrotron Radiation for MEdical Physics) beamline of the Italian synchrotron ELETTRA (Trieste, Italy).
Whole paraffin embedded samples were scanned with a sample to detector distance of 500 mm. A a 16-bit, water-cooled sCMOS camera (Hamamatsu C11440-22C ORCA-Flash 4.0 v2) was used to acquire 3600 angular distributed projections with an exposure time of 50 ms and an isotropic voxel size of 2.33 µm. Scans were performed in a 360° offset regime, resulting in a scan time of 180 s. A 0.5 mm Si filter was used, resulting in a mean beam energy of 19.6 keV. All SRµCT data sets were reconstructed using SYRMEP Tomo Project 48 software (STP Version 1.3.2). A single distance phase retrieval algorithm developed by Paganin et al. 49 was used with a delta over beta ratio of 100. To image the entire tumours, 3-4 single acquisitions were performed, and the resulting reconstructions were manually stitched together using Fiji (NIH).

Histology and immunohistochemistry (IHC).
Paraffin-embedded samples were sectioned in slices of 2 µm thickness using a standard microtome (HM 340 E microtome, Thermo Fischer Scientific). Histological haematoxylin and eosin (HE) and Masson trichrome staining (MTS), as well as immunostaining of SV40 T-antigen (T-Ag) marker and the macrophage specific marker CD68 were performed as described previously 29 . In brief, a Masson trichrome kit (MTS, Merck) was applied to stain collagen. A polyclonal rabbit anti-mouse CD68 antibody (abcam 125212, 1:500) was used to stain macrophages and monocytes. Since H8N8 breast tumour cells express T-Antigen, this protein was used as a marker for tumor cells and specifically stained with anti-SV40 T-Antigen antibody (homemade rabbit antibody, 1:5000) 50 .
Histological images were acquired using an Axiovert 200 inverted microscope (Zeiss). To depict the whole tumour sections, the built-in stitching algorithm of the microscope was used with a 10 × magnification.

Elastic registration of microCT and histological images.
Resampling of the CT dataset and identification of the appropriate matching virtual slice was done using VG Studio MAX (V3.1, Volume Graphics). Fiji was used to pre-process the histological images. Calculation of the deformation and transformation of the images was done via Elastix and Transformix 34,35 that were integrated into a tailored graphic user interface called Fuxlastix. The Fuxlastix frontend was developed in C++, using QT5.15 and is available on https:// github. com/ xPITc oding/ Fuxla stix. For the presented data, the registration process took less than a minute using a standard computer with a quad-core CPU.
Quantification of registration accuracy. The images acquired by CT and histology were compared in a block matching technique, calculating the mutual information MI and the displacement vector d for each block i. DI is then calculated as with G = i MI −1 i and d i the length of the displacement vector, as previously described 7 . Images with a high similarity are therefore characterised by a minimal DI. Fifteen tumour sections stained with histological staining www.nature.com/scientificreports/ protocols were evaluated. Statistical significance was determined using an unpaired, two-tailed Student's t-test. A p-value of 0.05 or lower was considered statistically significant.
Tissue classification using mixed density models. Classification of different structures in the combined data of histology and microCT was based on the assumption that after successful registration each pixel is characterised by four values-the three colour components of the histology image and the grey value of the microCT data set. Thus, similar tissue/cell types are characterised by similar combinations of those values. Classification was then performed by cluster analysis. Because the feature space is strongly overlapping due to partial volume effects and limitation of the registration approach, we used a mixed model of Gaussian distributions 51 . The parameters of these distributions were then evaluated by an expectation maximisation algorithm (EM) resulting in probabilities for each pixel to belong to a certain component/class of this model. We used the EM implementation of the python module scikit-learn 52 .