Xeno-learning: knowledge transfer across species in deep learning-based spectral image analysis

Data collection

The HSI animal data were collected at Heidelberg University Hospital following approval from the Committee on Animal Experimentation of the regional council Baden-Württemberg in Karlsruhe, Germany (approval nos. G-161/18, G-262/19 and G-62/23). The animals were treated in accordance with German laws for animal use and care, as well as the directives of the European Community Council (2010/63/EU).

The HSI human data were obtained during the SPACE trial (Spectral Characterization of Organs and Tissues during Surgery) at Heidelberg University Hospital, with approval from the Ethics Committee of the Medical Faculty of Heidelberg University, Germany (S-459/2020). The trial adhered to the ethical principles of the Declaration of Helsinki³⁰ and the principles of Good Clinical Practice³¹. The trial’s reporting followed the recommendations of the Consolidated Standards of Reporting Trials (CONSORT) guideline³². The SPACE trial was registered with the Research Registry (researchregistry6281) on 23 November 2020.

Hyperspectral image acquisition

The HSI camera system Tivita Tissue (Diaspective Vision GmbH) was used to collect the HSI data. This system captured hyperspectral images in a push-broom manner with a spectral resolution of approximately 5 nm, covering the spectral range from 500 nm to 1,000 nm. The resulting data cubes have dimensions of 640 × 480 × 100 (width in pixels × height in pixels × number of spectral channels). The camera system imaged an area of about 20 × 30 cm. An integrated distance calibration system, consisting of two light marks that overlap when the distance is correct, ensured an imaging distance of around 50 cm. The image acquisition process took approximately 7 s.

In addition to the HSI data cubes, the camera system computed functional parameters such as StO₂ (ref. ³). Furthermore, red–green–blue (RGB) images were reconstructed from the HSI data by combining spectral channels that capture red, green and blue light. More technical details on the hardware and the performed calculations can be found in cited publications³³.

To prevent spectral distortion from stray light, all other light sources were turned off during image capture and window blinds were closed. Motion artefacts were minimized by (1) mounting the camera on a swivel arm to keep it stationary during image capture, thus eliminating camera motion, and (2) capturing images from static scenes with no surgeon-induced object movements. Consequently, any motion artefacts would only be due to natural causes such as respiration and heartbeat. The camera perspectives were chosen to provide a clear view of all organs of interest in the scene.

For the pig and rat species, standardized recordings were carried out with a predefined image acquisition protocol¹¹. These data are used in Supplementary Fig. 3.

For the acquisition of animal data, the organs were mobilized and prepared to enable the spectral recording of representative tissue surfaces. Time series data were acquired by taking images every 30–60 s after the clamping procedure.

For the purpose of this work, the following types of visceral organ data were recorded: physiological (all three species), malperfused (all three species) and with contrast agent ICG (mice and porcine models). Note that that there was no opportunity to record ICG data in human patients for ethical reasons.

The animal malperfused organ data were acquired by vascular clamping of the infradiaphragmal aorta and caval vein to induce visceral malperfusion in a highly controlled and standardized manner. Haemodynamic stabilization and tissue oxygen desaturation over 2 min were awaited before measurements.

The human data were acquired in a clinical trial during which malperfused organs were recorded whenever they occurred intraoperatively, such as during kidney transplantation or acute organ ischaemia due to embolic or thrombotic events or iatrogenic preparation and resection. Human malperfused data were therefore available for the colon, liver, kidney and spleen. Consequently, analysis of malperfused data in the animal species was also restricted to this organ selection.

The animal ICG data were acquired by body weight-adapted intravenous application of commercial ICG. Pigs received 25 mg ICG injected into peripheral vein catheters, while rats received 2.5 mg directly injected into the caval vein. Pharmacological invasion and systematic distribution over 1 min were awaited before measurements.

Hyperspectral image annotation

From the 14,013 images, 2,596 were fully semantically annotated. For the remaining 11,417 images, polygon annotations were performed to annotate highly representative areas of the organ. All annotations were based on the reconstructed RGB image.

Polygon annotations were performed in the same manner as described previously², and the standardized recordings for the pig species are publicly available¹¹. Polygon annotations were used in the mixed effects analysis of Supplementary Fig. 3 and for several malperfused tissue annotations.

The semantic physiological porcine annotations are the same as described in ref. ⁸, where the semantic annotation process was performed by two different medical experts. Conflicts were revised by the same two medical experts. For all remaining semantic annotations, the annotation process was conducted by a team of medical experts using the Medical Imaging Interaction Toolkit³⁴. To ensure consistent labelling, all annotations were reviewed by the same medical experts. The segmentation networks were only trained on images with semantic annotations.

From the annotations, only the classes stomach, small bowel, colon, liver, pancreas, kidney, spleen, omentum, lung, skin, peritoneum and background were selected for this study, matching the organs which are available for all three species (see below). All remaining classes were not included in the analysis.

Data preprocessing

To mitigate sensor noise and transition the acquired HSI data from radiance to reflectance, the raw HSI data cubes were automatically calibrated using prerecorded white and dark calibration files by the camera system, as outlined in cited literature³. Calibration of the camera was performed before each surgery by taking a new white and dark image to compensate for various sources of signal distortion, such as attenuation effects of the light source³⁵.

After exporting the HSI cubes from the camera system, each pixel in the HSI cube was L1 normalized across the spectral channels to account for multiplicative illumination changes, such as those caused by fluctuations in the measurement distance.

Data statistics

The HSI database used in this study comprises 14,013 images from 319 subjects in three species. Annotation has been performed for 12 classes (11 organs and background). These organs were selected as they were available for all three species owing to anatomical, anaesthesiologists and technical considerations. Malperfused tissues were captured on 446 pig images, 1,538 rat images and 796 human images. Tissues following ICG injection were captured on 564 pig images and 1,137 rat images. The remaining 4,019 pig images, 4,683 rat images and 830 human images show physiological tissues. A detailed overview of the database is shown in Supplementary Fig. 1.

Segmentation models

The used segmentation networks are the same as described previously⁸, with the organ transplantation extension proposed^9,10. In short, the full HSI data cube is passed on to a U-Net with an efficientnet-b5 encoder pretrained on the ImageNet dataset. Dice loss and cross-entropy loss are equally weighted and computed for all valid pixels inside a batch, that is, every pixel which does not belong to one of the ignored classes.

The same hyperparameters were used for all models. Adam³⁶ was used as an optimization algorithm with an exponential learning rate scheme (initial learning rate η = 0.001, decay rate, γ = 0.99, Adam decay rates β₁ = 0.9 and β₂ = 0.999). Training was carried out with a batch size of 8 for 100 epochs, with each epoch consisting of 500 images. During the last ten epochs, stochastic weight averaging was applied³⁷. Underrepresented classes were oversampled to ensure an equal class distribution.

During training, images were augmented to increase the size and the diversity of the training data. First, the same affine transformations as described previously⁸ (shift, scale, rotate and flip operations) were applied. Then, target tissues in the images (kidney in our study) were transformed with our proposed xeno-learning method. Finally, the organ transplantation method proposed in refs. ^9,10 followed by a L1 re-normalization of the image was applied.

Training and validation set-up

A similar training and validation set-up was used for all networks. Validation was carried out on the basis of a nested cross-validation scheme to provide a more robust performance estimation based on the entire dataset³⁸. The number of outer folds was set to 3 and the number of inner folds was set to 5. The folds were generated on the basis of iterative stratification for multilabelled data³⁹ to ensure a similar label distribution across folds. Final predictions for an image were obtained by ensembling the softmax output from all available networks.

To prevent model overfitting in the standard ML approach, it is crucial to evaluate methods using an untouched test set. We extend this principle by incorporating completely unseen tasks and species in our evaluation framework. During the development of our method, we focused exclusively on one task (malperfusion) and data from only the pig species, phrasing the problem as a pig-to-pig malperfusion task, where the domain shift stems from unseen individuals. After finalizing our method, we first transferred the approach to further species (rats and humans) and then to another knowledge-transfer task. Crucially, pig ICG data, as well as all rat and human data, were not used at any stage during the development of our method.

Following the recommendations from previous publications⁴⁰ and to overcome the limitations of individual metrics, we assessed the segmentation performance via the overlap-based DSC and the boundary-based NSD. The class-specific thresholds for the NSD were set to the same values as reported previously⁸. Metric scores were always calculated per class and then hierarchically aggregated towards a final class-level score. That is, class scores were averaged first across all images of the same subject and then across all subject-level class scores.

For the results of Fig. 6, only the DSC is reported because some of the annotations of malperfused organs were performed via polygon annotations (see ‘Hyperspectral image annotation’ section). The goal of these polygon annotations was to annotate representative areas of an organ but not to strictly follow the organ boundaries. Hence, the calculation of the NSD (which assesses boundary agreement) was not reasonable.

Confidence intervals reported in this study reflect the subject-level sampling variability of the data. After the aggregation towards subject-level scores, bootstrapped sampling was performed with replacement 1,000 times for each organ, and the sampled subject-level scores were averaged to retrieve an organ-level score per bootstrap sample. After aggregating (for example, mean or median) across organs for each bootstrap sample, 95% confidence intervals were calculated.

Linear mixed model analysis

Separate linear mixed models for each wavelength and organ were employed to analyse explained variation to evaluate the relevance of factors contributing to changes in the observed spectrum (Supplementary Fig. 3). The proportion of explained variance was derived through the empirical decomposition of explained variation on the basis of the variance components version of the mixed model⁴¹.

More precisely, linear mixed models were fitted for each organ and wavelength separately, with fixed effects for the factor angle and the factor species as well as random effects for the factor subject and the factor image.

$${\mathrm{reflectance}}_{i{jk}}=\alpha +{\mathrm{species}}_{i{jk}}^{\top }{\rm{\cdot }}\beta +\,{\mathrm{angle}}_{i{jk}}^{\top }{\rm{\cdot }}{\boldsymbol{\theta }}+{\delta }_{i}+{\gamma }_{ij}+{\varepsilon }_{i{jk}},$$

(1)

for repetition k = 1,…, 3 of image j = 1,…, n_i of animal i = 1,…, 24 (11 pigs and 13 rats).

The number of images n_i varied per animal and organ. Here, α denotes a fixed intercept, ${\mathrm{species}}_{i{jk}}^{\top }$ is a row vector of length 2 indicating the species rat or pig and β denotes the corresponding fixed effect. Similarly, θ is a vector of fixed effects corresponding to the camera angles (‘perpendicular to tissue surface’, ‘25 degrees from one side’ and ‘25 degrees from the opposite side’). The random intercept δ_i ~ N ${(0,\,{\sigma }_{\delta }^{2})}_{\gamma }$ describes animal-specific variation, and the random intercept γ_ij ~ N ($0,\,{\sigma }_{\gamma }^{2}$) describes image-specific variations. The residuals ε_ijk ~ N ($0,\,{\sigma }_{\varepsilon }^{2}$) capture the variability between repeated recordings of the same image. Within the model, we assumed that the random effects and the residuals are stochastically independent.

In addition, 95% pointwise confidence intervals were obtained on the basis of parametric bootstrapping with 500 replications for an indication of the uncertainty in the relevance estimates.

Xeno-learning

The driving insight of our concept is that shared (for example pathophysiological) mechanisms manifest in comparable relative spectral changes. In our approach (Fig. 2), we leverage this insight to transfer knowledge across species through data augmentation. Specifically, we first learn the effect of certain interventions, such as clamping or contrast agent injection, in the source species and then apply this knowledge in the target species.

Learning relative changes in the source species

We encode relative spectral changes in a transformation matrix that can be applied to any species. For the specific case of perfusion, a set of linear transformations t_i(s_p) are learnt that transform physiological spectra s_p to malperfused spectra s_m (the transformation is applied independently for each spectrum). To cover a variety of different perfusion states, we learn a whole set of transformations, each of which represents the spectral change between physiological and malperfused kidneys. Each transformation is a linear model represented by two parameters: a weight matrix W_i $\in$ R^{100 × 100} and a bias vector b_i $\in$ R¹⁰⁰, by analogy with a multivariate linear regression model, so that the transformation is defined as

$${t}_{i}({s}_{{\rm{p}}})={W}_{i}{\rm{\cdot }}{s}_{{\rm{p}}}+{{\bf{b}}}_{i}.$$

(2)

To learn the transformation t_i, image pairs (p_i, m_i) are randomly selected in the source species, consisting of an image p_i showing a physiological kidney and an image m_i showing a malperfused kidney. The images do not necessarily need to come from the same subject (the transformation will always be applied to different subjects in the target species).

The transformation parameters W_i and b_i should be based on the spectra in the image pair (p_i, m_i). However, because the number of spectra for a specific organ in each image is usually different (and hence no match between the spectra s_p and s_m is available), W_i and b_i cannot be optimized directly in closed form. Instead, the following indirect optimization scheme is applied for each selected image pair (p_i, m_i):

(1)

The parameter weight matrix is initialized with the identity matrix W_i = I and the bias with zeros b_i = 0. Hence, the initial transformation t_i(s_p) does not modify the spectra.
(2)

The MSE loss is used to compare the malperfused spectra s_m, $\in$ S_m (set of all spectra in the malperfused image m_i) with the physiologically transformed spectra t_i(s_p) with s_p $\in$ S_p (set of all spectra in the physiological image p_i). The loss consists of three components: (1) two histograms each with 50 bins; h_m calculated from the normalized reflectance values from all s_m and h_t calculated from the normalized reflectance values from all t_i(s_p), (2) the mean spectrum sˆ_m and tˆi(s_p) and (3) the s.d. spectrum s˜_m and t˜_i(s_p).
(3)

The linear model is fitted iteratively in 100 steps with the Adam optimizer (learning rate η = 0.001, Adam decay rates β₁ = 0.9 and β₂ = 0.999) while gradually adapting the parameters W_i and b_i.

Transferring knowledge

After the optimization phase in the source species, all learnt changes, encoded in the parameters W_i and b_i can be applied in the target species. For increased data variety, this is done dynamically during the training process as augmentation. For every image in the batch, the physiological kidney spectra S_p are selected using the available segmentation mask. Then, a parameter set j is randomly selected and the corresponding parameters W_j and b_j are used to transform every kidney spectrum s_p $\in$ S_p using equation (2). To cover many different perfusion states, the transformed spectra are linearly interpolated with a randomly selected weight λ $\in$ (0; 1)

$$s=(1-\lambda ){\rm{\cdot }}{s}_{{\rm{p}}}+\lambda {\rm{\cdot }}{t}_{{\rm{j}}}({s}_{{\rm{p}}}),$$

(3)

where s denotes the replaced spectra which are shown to the network. In total, this augmentation is applied to an image with probability P = 0.8 so that some unaltered kidney spectra are also shown during the training process.

Distribution comparison

The extended training distribution shown in Fig. 7 was derived by a training process comprising 100 epochs with 500 images per epoch, during which our data augmentation was applied in the same manner as during the training of the segmentation networks. For each organ in each image seen during the training phase, the median spectrum was computed before and after augmentation. The resulting distributions of median spectra are referred to as baseline training distribution (Fig. 7, black distribution) and extended training distribution (Fig. 7, purple distribution).

To compare the baseline training distribution with the extended training distribution, two-dimensional projections of all median spectra were computed using PCA, where the principal components were derived from the baseline training distribution. The median spectra from the extended training distribution were subsequently projected into this same feature space. Both distributions were modelled using kernel density estimation. Distances were calculated by measuring the Euclidean distance between each median spectrum in the test set and the nearest median spectrum in either the baseline or extended training distribution. These distances were then hierarchically averaged, first across all images belonging to a single subject and subsequently across all subjects.