DREAMER-S: Deep leaRning-Enabled Attention-based Multiple-instance approaches with Explainable Representations for Spatial biology.

Multiple Instance Learning (MIL) is an innovative machine learning approach designed to address challenges associated with weakly labelled data, where class labels are assigned to collections of instances, known as bags, rather than to individual instances themselves. A bag is labelled positive if at least one of its instances is positive, while it is labelled negative only if all instances are negative [ 1 ]. This approach is particularly beneficial in scenarios where labelling is expensive or impractical, allowing for the classification of bags based on the presence of positive instances within them [ 2 , 3 ]. Hence, it is a common strategy in in digital pathology, to analyse whole slide images (WSIs) which are subdivided into smaller patches so that a deep neural network model can be trained for specific purposes such as detecting epithelial cells, grading tumours, and identifying other histopathological features (reviewed systemically in [ 4 ]). MIL can therefore facilitate the simultaneous feature localisation and classification of disease in a computationally-efficient manner in large images, addressing the inherent challenges posed by the complexity and size of WSIs [ 5 ]. Mid-infrared (IR) spectroscopy allows for the identification of specific molecular signatures associated with various pathological conditions using comprehensive biochemical profiles of tissues through the analysis of the absorption of infrared light. This is because biomolecules absorb in the mid-IR (4000 cm −1 to 400 cm -1 ) producing an absorbance spectrum which can be treated as a unique biochemical fingerprint. This fingerprint can be used to differentiate between normal and tumour tissues [ 6 ], identify spectral signatures indicative of endometriosis [ 7 ], elucidate the chemical composition of lipid droplets by size in liver tissue and endothelial cells [ 8 ], amongst other applications. IR spectroscopy is therefore an attractive approach for histopathological analysis due to its non-invasive, label-free, objective, and quantitative nature enabling clinically informed decisions without the need for extensive sample preparation or labelling [ 9 ]. In a previous study, from which the present study stems, [ 10 ], the response of two patient derived xenograft (PDX) mouse colon cancer models were studied (CRC0344 and CRC0076), which differed in their responsiveness to chemotherapy. The CRC0344 model was found to be sensitive to 5-fluorouracil (5-FU)-based chemotherapy, while the CRC0076 model was less responsive to chemotherapy but responded with tumour regression by addition of the apoptosis sensitiser and BCL-2 antagonist, ABT-199. Traditional imaging methods such as 18F-Fluorodeoxyglucose positron emission tomography/computed tomography (18F-FDG-PET/CT) were demonstrated to have potential as early response biomarkers [ 10 ]. The central hypothesis of the present study was that applying a deep learning MIL (DL-MIL) approach to the rich biochemical data from IR chemical imaging of the PDX colon cancer tissue would uncover unique spectral biomarkers capable of discriminating each PDX-model on their sensitivity to chemotherapeutic treatment. Modern IR imaging microscope systems employing focal plane array (FPA) detectors, or similar, allow the generation of full spectral information spatially across a sample, such that the output “chemical” images are termed hyperspectral data cubes. The challenge with using these images within a deep learning modelling approach in the absence of dimensionality reduction is the sheer size of each image in data terms. For example, a single chemical image generated by a 256 × 256 FPA array capturing 217 wavenumbers per spectrum will contain 14,221,312 datapoints. Comparatively a typical image size within the CIFAR10 dataset commonly used as a training dataset for deep learning only contains 32 x 32 spatial pixels, each with 3 RGB channels, or 3072 datapoints [ 11 ]. Therefore, our hypothesis is that the use of a deep learning-MIL framework could represent a practical approach in this context, in terms of both computational efficiency and identification of spectral region of interests (ROIs) with biochemistry that are linked to the classification target (in this case identification of chemotherapeutic resistance or sensitivity). A number of recent studies have employed weakly supervised approaches with vibrational spectroscopic data for the analysis of both imaging or point-spectral data. For example, MIL and Support Vector Machine (SVM) algorithms have been applied to Raman spectroscopy for the detection of COVID-19 in saliva samples [ 12 ], while MIL and Random Forest (RF) algorithms have been utilised in IR studies to identify aberrant tissue in mouse livers [ 13 ]. Research by Shi et al [ 14 ] and Phan et al [ 15 ] has demonstrated the potential of the MIL approach with CNNs for the identification of regions of interest with spectral images containing signatures of drug fingerprints and microplastics, respectively. While MIL-DL frameworks have been extensively adopted in other biological domains like digital pathology [ 16 ] and transcriptomics [ 17 ], their potential for application to spatial biological imaging datasets remains largely untapped [ 18 ]. To date, there is a lack of focus on the development of an integrated MIL-DL architecture incorporating explainability techniques to pinpoint the specific spatial biomarkers that underpin systemic responses to treatment in heterogeneous tissue environments. The primary objectives of this study are to establish a weakly supervised MIL-DL framework for this purpose, with an exemplar application to QCL-IR chemical imaging data in the context of the classification of chemotherapeutic sensitivity within a PDX-model, and to determine the optimal model architecture for this application. To achieve this, a two-step grid-parameter search was utilised to tune hyperparameters, including network depth, hidden nodes, residual connections, and learning rates. Furthermore, the study evaluates the decision-making processes of the best-performing model by employing Explainable AI (XAI) - specifically attention mechanisms and SHAP values - to identify and extract candidate spectral biomarkers of chemotherapeutic sensitivity. These features are further explored via UMAP clustering to correlate learned spectral signatures with biological relevance, specifically protein expression. Ultimately, while demonstrated here in spectral histopathology, the DREAMER-S framework serves as a versatile data-reduction and ROI-filtering tool for high-content spatial biology pipelines such as Raman/CARS and spatial transcriptomics where manual pixel-level annotation is typically infeasible.

This study successfully demonstrates the development and application of an integrated deep learning-multiple instance learning (DL-MIL) framework for downstream classification tasks with high-content imaging data. In the present instance we apply this approach to the classification of chemotherapy sensitivity using QCL-IR hyperspectral imaging of colorectal cancer PDX models. Through a systematic optimisation process, we identified a robust, non-residual linear neural network architecture capable of distinguishing between chemo-sensitive (CRC0344) and chemo-resistant (CRC0076) models with high fidelity (F1-score = 0.95). The MIL attention mechanism proved highly effective, autonomously identifying salient tissue regions from over 200,000 spectral instances per image, thereby circumventing the need for laborious manual annotation. Explainable AI, through the use of SHAP, provided unprecedented insight into the model’s decision-making process as it reveals that the classification was not based on arbitrary patterns but on distinct and biochemically relevant spectral features, primarily related to nucleic acid phosphate backbones, lipid carbonyl groups, and protein amide bands. Critically, this analysis also demonstrated the model’s sophistication in learning to ignore confounding signals from experimental artifacts such as paraffin wax. Furthermore, we validated the biological relevance of the model-attended spectra through two distinct approaches. Unsupervised UMAP clustering of these spectra successfully stratified the samples not only by their patient-of-origin but also by their specific treatment group, confirming that the learned features were rich with pharmacodynamic information. Secondly, correlation studies established a direct, albeit highly context-dependent, link between these spectral markers and the expression of key apoptotic proteins, Bim and Puma. The discovery of distinct biomarker pairs for each PDX model (e.g., Amide I/Puma for CRC0076 vs. Amide III/Puma for CRC0344) underscores the importance of discovering patient-specific molecular signatures. In summary, this research presents a powerful, end-to-end pipeline for the label-free, data-driven discovery of spatial biomarkers from complex high-content imaging data by leveraging the entire spatial signature without prior feature engineering. It also highlights the potential of explainable AI-driven spectral histopathology to advance pre-clinical drug efficacy studies and paves the way for developing more precise, patient-specific diagnostic and prognostic spectral-driven tools for biomedical and pre-clinical applications in the future.

Here, we utilised CRC PDX material previously generated by O’Farrell et al. during a preclinical study [ 10 ], which was approved by the Health Products Regulatory Authority [HPRA; #AE18982-P099] and the University College Dublin Animal Research Ethics Committee [AREC; #AREC-16–11]. PDX material for that study was collaboratively provided by researchers at the Candiolo Institute, Torino group from the PDX bioresource [ 19 , 20 ], with coded/de-identified material shared. Liver metastasectomy samples used to establish the original PDX models were obtained under informed consent from patients at the Candiolo Cancer Institute (Candiolo, Torino, Italy), Ospedale Mauriziano Umberto I (Torino), and Città della Salute e della Scienza di Torino – Presidio Molinette (Torino). Samples were procured under the approval of the Review Boards of the Institutions (PROFILING protocol No. 001-IRCC-00IIS-10). Two PDX models were selected based on their predicted sensitivity to chemotherapy, one sensitive to 5-FU-based chemotherapy (labelled as CRC0344) and another sensitised by ABT-199 (labelled as CRC0076). Animals bearing either subcutaneous CRC0076 or CRC0344 tumours received one of the following treatment regimens: (i) FOLFOX (all constituent drugs delivered IP, once weekly (on day 3 of each cycle), 5-FU (40 mg/kg in PBS) + Folinic Acid (13.4 mg/kg in PBS) followed two hours later by OX (2.4 mg/kg in 5% glucose/water (v/v)), (ii) ABT-199 (oral gavage, once daily at 100 mg/kg, dissolved in 60% phosal, 50 propylene glycol (PG), 30% polyethylene glycol (PEG) 400 and 10% ethanol), (iii) FOLFOX with ABT-199 (dosed as described above) and (iv) a vehicle control (all diluent vehicle solutions), resulting in eight treatment classes. Researchers undertaking the animal studies were blinded to the expected sensitivity of the PDX models. At the end of the study (i.e., after 4 weeks of treatment or when the study reached humane end point [based on a scoring system accounting for tumour size and animal wellbeing]), tumours were excised, rinsed twice in Dulbecco’s-(D)PBS and fixed in 4% formaldehyde for 48 h and embedded in paraffin. A 5 µm-thick tissue section was prepared and mounted on a calcium fluoride (CaF 2 ) slide for analysis. To retain all tissue the section was not subjected to any chemical dewaxation. Spectroscopic measurements were carried out using a Daylight Spero-QT 340 Quantum Cascade Laser (QCL) infrared microscope operating in transmission mode. Hyperspectral chemical images (HCIs) were captured at a low magnification (0.3 NA) across the wavenumber range 952–1800 cm −1 . The generated HCIs comprised a spatial area of 480 × 480 pixels and a spectral depth of 213 wavenumbers. Following data acquisition, a multi-step pre-processing pipeline was applied to the raw spectral data using custom scripts written in Python (version 3.10). The primary objective of these steps was to correct for artefacts and to isolate the spectra originating specifically from the tissue. First, a rubber-band baseline correction was applied to each spectrum using the Pybaselines library (version 1.2) to account for scattering effects [ 21 ]. To distinguish viable tissue regions from the slide background, a fixed absorbance threshold of 0.1 was applied at the Amide I peak (located at approximately 1654 cm -1 ). This threshold was determined empirically to ensure the exclusion of low-signal background noise while preserving tissue pixels. Spectra failing to meet this criterion were classified as non-tissue and replaced with zero-vectors to maintain the spatial dimensions of the image. Finally, to account for variations in tissue thickness, vector normalisation was performed on the remaining tissue-positive spectra using Scikit-Learn module (version 1.7.0) [ 22 ]. The full dataset for this study comprised 40 hyperspectral images. Each image, containing 213 wavenumbers and 230,400 individual spectra, represents over 49 million data points in each image. The complete pre-processing workflow, from raw data to corrected spectra, is visually summarised in Fig 1a . (a) A summary of the preprocessing steps for chemical images, including preparation procedures prior to training the deep learning model and optimising hyperparameters. (b) The skeletal structure of the linear model architecture utilised in this study, which is loosely inspired by the ResNet framework. For each chemotherapy sensitivity group (CRC0076 and CRC0344), the hyperspectral image dataset was manually partitioned into equal subsets (50%:50%) at the animal (mouse) level rather than the image level for training and validation. Each subset comprised of 20 hyperspectral images: 10 hyperspectral images for each of CRC0076 and CRC0344 ( S1 Table ). Both training and validation datasets contained a representative distribution of the treatment conditions described earlier (FOLFOX, ABT-199, FOLFOX with ABT-199 and a vehicle control) and each group was separated by the identity of individual mice to prevent data leakage during model training and testing. In our multiple instance learning (MIL) approach, each hyperspectral image file was treated as a bag consisting of multiple spectral instances (230,400 individual spectra in each bag). The  .MAT files were converted into spectral tensors that has shape of (Batch number, Spectral instances, Features); (1, 230400, 213) in a custom Pytorch dataset class. Each bag was assigned a label corresponding to its respective chemotherapy sensitivity group (label 0 for CRC0076 and label 1 for CRC0344) for the classification model training task. During the MIL training, the whole bag (chemical image) was fed through the network for a single prediction class. Deep learning model architectures were developed using Pytorch (version 2.4.0) [ 23 ] with their design loosely inspired by the Residual-Net (ResNet) block where each block consists of repeating linear, batch normalisation and ReLU activation layers (see Fig 1b ) [ 24 ]. An optional residual connection was introduced before the first linear layer that connected the final ReLU activation layer in each block. An input layer was also introduced prior to the first block to allow the input of spectral data, and a final linear layer was introduced for multi-class classification. Following the sequence of blocks, and preceding the final classification layer, a MIL-attention mechanism was applied based on the implementation from [ 25 ]. The purpose of this mechanism is to produce a single, representative embedding for the entire input “bag” by weighting the significance of each instance. This allows the model to identify which instances are most influential for the final prediction. Subsequently, a final linear layer uses this attention-pooled embedding for multi-class classification. The MIL attention mechanism can be summarized as follows: let H={h 1 , h 2 , …, h K } represent a bag containing K instance embeddings, where each instance h k  ∈ R M is an M -dimensional real-valued feature vector representing an individual instance. compute a weighted average, z , of these instances to form a single bag-level representation: The attention weights, a k , are calculated using a small neural network and normalized across all instances in the bag via the softmax function. This ensures that the weights sum to 1 and allows the model to learn the relative importance of each instance: where V  ∈  R L×M and w  ∈  RL are the learnable weight parameters of the attention network. This formulation enables the model to focus selectively on the most informative instances within each bag, thereby enhancing the overall discriminative power of the bag-level representation in weakly supervised learning settings. To determine the optimal model configuration for our hyperspectral data, we conducted a systematic grid search of key hyperparameters. This optimisation process focused on two main areas: 1) the model’s fundamental architecture and 2) its learning rate parameter. The main goal was to maximise classification performance while balancing model interpretability, maintaining computational efficiency, and mitigating the risk of overfitting. The architectural search explored several core parameters. First, network depth was varied by adjusting the number of sequential blocks (N B ), evaluating shallow (N B  = 1) and deeper (N B  = 2 or 3) networks. Second, the network width was modified by changing the number of nodes (N HLN ) within the fully connected layers, with values ranging from 64 to 256. We also investigate the network’s expansion strategy by assessing two distinct approaches. In the ‘expanding’ configuration, the number of nodes was doubled in each successive block, a design intended to capture increasingly complex and hierarchical features. Conversely, the ‘shrinking’ configuration halved the number of nodes per block, progressively reducing the network’s complexity. Furthermore, the inclusion of residual connections was tested as an additional parameter to assess their potential to improve gradient flow and overall model performance, particularly in deeper architectures. The optimisation was performed in a two-stage process. Initially, the optimal architectural parameters were identified using a default learning rate (1x10 -3 ). Once the best-performing architecture was established, its specific learning rate for the Adam optimiser was then fine-tuned. Throughout this process, Cross-Entropy was employed as the loss function, while the F1-score served as the primary metric for evaluating model performance during the training and validation phases. To ensure the robustness of our findings, each hyperparameter combination was trained for 20 epochs and the process was repeated at least three times using different randomisation seeds within the Pytorch module. The final selection was based on the parameter set that yielded the highest average F1-score on the validation dataset. A comprehensive summary of the hyperparameters and their tested values is provided in Table 1 . Each hyperparameter set was evaluated across at least three independent runs, and the configuration with the highest F1-score on the validation set was selected. This chosen configuration was then retrained for 50 epochs, and the best-performing model checkpoint was saved for further analysis. To gain insight into the model’s decision-making process, we analysed the attention weights assigned to the spectra within the validation dataset. This attention mechanism is a key feature for interpretability, as it reveals which spectral features and spatial regions the optimised model considered most significant for classification. First, the validation dataset was processed by the fully trained model, and the attention weight for each spectral instance was extracted. For each hyperspectral image, these weights were aggregated and then normalised to a range of 0–1. To visualise the spatial distribution of these important regions, the normalised attention weights were reshaped into two-dimensional ‘attention maps’ matching the original chemical image dimensions (480 x 480 pixels). These maps were then overlaid onto a greyscale image representing the tissue morphology region on the slide. This morphological base image was generated by calculating the Area Under the Curve (AUC) for each spectrum, which provides the total integrated signal intensity at each pixel ( i,j ). The AUC was calculated using the following integral, where P (ν) is the spectral intensity at wavenumber, ν : Finally, for more targeted downstream analysis, we identified and isolated the most influential spectra as determined by the model. From each chemical image in the validation set, the spectra corresponding to the top 5% of attention weights were filtered and aggregated. This process created a refined dataset comprising only the spectra that the model deemed most salient, providing a high-value subset for subsequent investigation. To gain a deeper insight into the spectral features driving the model’s predictions, we analysed the data using SHapley Additive exPlanations (SHAP). This analysis was performed with the shap Python library (version 0.48.0) [ 26 ]. The primary objective was to identify the key wavenumbers that have the most significant impact on distinguishing between classes, thereby enhancing model interpretability and facilitating chemical inference from its outcomes. For this task, we specifically employed SHAP’s GradientExplainer, an approach optimised for deep learning models. This explainer was applied to the validation dataset to calculate feature importance scores. The resulting SHAP values were then used to generate a summary plot for each class ranking the wavenumbers based on their overall importance to the model and illustrate not only the magnitude, but also the direction of each feature’s effect. This will show whether a high absorbance at a specific wavenumber pushes the prediction towards or away from a particular class. It also allows for the direct identification of the key spectral bands that underpin the model’s performance. To determine if the top spectra identified by the model’s attention mechanism possess biological relevance, we investigated whether they could be separated based on treatment type. This analysis was performed using Uniform Manifold Approximation and Projection (UMAP), implemented via the umap-learn Python library (version 0.5.7) [ 27 ]. UMAP is a powerful non-linear dimensionality reduction technique used to project high-dimensional data into a low-dimensional embedding, making it easier to visualise data structure and class separability [ 27 ]. Prior to UMAP analysis, these input spectral features were scaled using a Standard scaler from the scikit-learn library [ 20 ]. This pre-processing step normalises the distribution of each feature, ensuring that the intensity scaling is consistent across all wavenumbers and improving the robustness of the subsequent distance-based calculations in UMAP. A systematic parameter search was conducted to find the optimal UMAP configuration that could most effectively cluster the spectra according to their treatment class (FOLFOX, ABT-199, FOLFOX with ABT-199, or vehicle control). The key parameters explored included the number of neighbours ( n_neighbors , ranging from 5 to 80), the minimum distance between embedded points ( min_dist , ranging from 0.1 to 0.5), and various distance metrics (Euclidean, Cosine, and Manhattan). For each set of parameters, a two-dimensional UMAP embedding of the spectra was generated. This low-dimensional representation was then clustered using the k-means algorithm, with the number of clusters ( k ) set to seven, corresponding to the number of distinct treatment labels. The performance of each UMAP configuration was quantified by evaluating the agreement between the k-means cluster assignments and the ground-truth treatment labels. This was measured using the Normalised Mutual Information (NMI) score, a metric that assesses the similarity between the clustering and the ground truth class membership [ 28 , 29 ]. The NMI between the true labels ( U ) and the k-means cluster assignments ( V ) is calculated as: where, I(U,V) represents the mutual information between the two assignments, while H(U) and H(V) are their respective entropies. The optimal UMAP parameters were selected based on the configuration that yielded the highest NMI score, as this indicates the clearest and most meaningful separation of the treatment groups in the low-dimensional space. An additional correlation analysis was conducted to investigate the relationship between spectral data and biological endpoints. We aimed to investigate whether the changes of the pro-apoptotic proteins Puma and Bim, key mediators of chemotherapy-induced apoptosis and evaluated in the original study [ 10 ], lead to detectable spectral bands, using selected wavenumber and ratios. We extracted four protein-related spectral bands: the integrated intensity of the Amide I (1600–1700 cm −1 ), Amide II (1540–1560 cm  −1 ), and Amide III (1200–1350 cm  −1 ) bands, and the Amide I/II peak ratio (~1658/1544 cm  −1 ). Additionally, a series of spectral band ratios were also calculated to assess the relative content of proteins to other key biomolecules, including lipids and nucleic acids. These included ratios of Amide I to Amide II, DNA, RNA, and lipids, providing further insight into the overall biochemical changes within the cells. Details on the ratio’s used and their wavenumber bands are available in S2 Table . Then, using relative protein expression data acquired in the original PDX study [ 10 ], a Pearson correlation analysis was conducted independently for each chemotherapy sensitivity group (CRC0076 and CRC0344) to link the spectral metrics with apoptotic protein levels using correlation function available in Pandas python module (version 2.2.2).

(DOCX) (DOCX) (DOCX)