Leo Fillioux

Joseph Boyd

Maria Vakalopoulou

Paul-Henry Cournède

Stergios Christodoulidis

Structured State Space Models for Multiple Instance Learning in Digital Pathology

Keywords: Multiple instance learning, Whole slide images, State space models

Multiple instance learning is an ideal mode of analysis for histopathology data, where vast whole slide images are typically annotated with a single global label. In such cases, a whole slide image is modelled as a collection of tissue patches to be aggregated and classified. Common models for performing this classification include recurrent neural networks and transformers. Although powerful compression algorithms, such as deep pre-trained neural networks, are used to reduce the dimensionality of each patch, the sequences arising from whole slide images remain excessively long, routinely containing tens of thousands of patches. Structured state space models are an emerging alternative for sequence modelling, specifically designed for the efficient modelling of long sequences. These models invoke an optimal projection of an input sequence into memory units that compress the entire sequence. In this paper, we propose the use of state space models as a multiple instance learner to a variety of problems in digital pathology. Across experiments in metastasis detection, cancer subtyping, mutation classification, and multitask learning, we demonstrate the competitiveness of this new class of models with existing state of the art approaches. Our code is available at https://github.com/MICS-Lab/s4 digital pathology.

Introduction

Precision medicine efforts are shifting cancer care standards by providing novel personalised treatment plans with promising outcomes. Patient selection for such treatment regimes is based principally on the assessment of tissue biopsies and the characterisation of the tumor microenvironment. This is typically performed by experienced pathologists, who closely inspect chemically stained histopathological whole slide images (WSIs). Increasingly, clinical centers are investing in the digitisation of such tissue slides to enable both automatic processing as well as research studies to elucidate the underlying biological processes of cancer. The resulting images are of gigapixel size, rendering their computational analysis ⋆ These authors contributed equally to this work. challenging. To deal with this issue, multiple instance learning (MIL) schemes based on weakly supervised training are used for WSI classification tasks. In such schemes, the WSI is typically divided into a grid of patches, with general purpose features derived from pretrained ImageNet [START_REF] Russakovsky | Imagenet large scale visual recognition challenge[END_REF] networks extracted for each patch. These representations are subsequently pooled together using different aggregation functions and attention-based operators for a final slide-level prediction.

State space models are designed to efficiently model long sequences, such as the sequences of patches that arise in WSI MIL. In this paper, we present the first use of state space models for WSI MIL. Extensive experiments on three publicly available datasets show the potential of such models for the processing of gigapixel-sized images, under both weakly and multi-task schemes. Moreover, comparisons with other commonly used MIL schemes highlight their robust performance, while we demonstrate empirically the superiority of state space models in processing the longest of WSI sequences with respect to commonly used MIL methods.

Related work

Using pretrained networks for patch-wise feature extraction is a well established strategy for histopathology analysis [START_REF] Coudray | Classification and mutation prediction from non-small cell lung cancer histopathology images using deep learning[END_REF][START_REF] Tellez | Neural image compression for gigapixel histopathology image analysis[END_REF]. An extension of this approach is with MIL, where the patch-wise features of an entire slide are digested simultaneously by an aggregator model, such as attention-based models CLAM [START_REF] Lu | Data-efficient and weakly supervised computational pathology on whole-slide images[END_REF] and TransMIL [START_REF] Shao | Transmil: Transformer based correlated multiple instance learning for whole slide image classification[END_REF], the latter being a variant of self-attention transformers [START_REF] Vaswani | Attention is all you need[END_REF]. [START_REF] Chen | Scaling vision transformers to gigapixel images via hierarchical self-supervised learning[END_REF] proposes another transformer-based method in the form of a hierarchical ViT. Similar to our multitask experiments, [START_REF] Gao | A semi-supervised multi-task learning framework for cancer classification with weak annotation in whole-slide images[END_REF] explores combining slide-level and tilelevel annotations with a minimal point-based annotation strategy. One of the key components of MIL methods is the aggregation module that pools together the set of patch representations. Mean or max pooling operations are among the simplest and most effective for aggregating predictions over a whole slide [START_REF] Campanella | Clinical-grade computational pathology using weakly supervised deep learning on whole slide images[END_REF]. In contrast, recurrent neural networks (RNN) with long short-term memory (LSTM) [START_REF] Hochreiter | Long short-term memory[END_REF] model the patches more explicitly as a set of tokens in sequence. In particular, LSTM networks have been shown to work well in different MIL settings including both visual cognition [START_REF] Wang | In defense of lstms for addressing multiple instance learning problems[END_REF] and computational pathology [START_REF] Agarwalla | Representation-aggregation networks for segmentation of multi-gigapixel histology images[END_REF].

The state space model is a linear differential equation,

ẋ(t) = Ax(t) + Bu(t) y(t) = Cx(t) + Du(t) (1)
that is widely studied in control theory, and describes a continuous time process for input and output signals u(t) ∈ R p and y(t) ∈ R q , and state signal x(t) ∈ R n , and where the process is governed by matrices A ∈ R n×n , B ∈ R n×p , C ∈ R q×n , D ∈ R q×p . In HiPPO [START_REF] Gu | Hippo: Recurrent memory with optimal polynomial projections[END_REF] (high-order polynomial projection operator), continuous time memorisation is posed as a problem of function approximation in a Hilbert space defined by a probability measure µ. For a scaled Legendre probability measure, one obtains the HiPPO matrix A, which enforces uniform weight in the memorisation of all previously observed inputs, in contrast to the exponentially decaying weighting of the constant error carousel of LSTMs [START_REF] Hochreiter | Long short-term memory[END_REF]. The HiPPO mode of memorisation is shown empirically to be better suited to modeling long-range dependencies (LRD) than other neural memory layers, for which it serves as a drop-in replacement.

Whereas in HiPPO, the state matrix A is a fixed constant, the linear state space layer (LSSL) [START_REF] Gu | Combining recurrent, convolutional, and continuous-time models with linear state space layers[END_REF] incorporates A as a learnable parameter. However, this increased expressiveness introduces intractable powers of A. In [START_REF] Gu | Efficiently modeling long sequences with structured state spaces[END_REF], the LSSL is instead reparameterised as the sum of diagonal and low-rank matrices, allowing for the efficient computation of the layer kernel in Fourier space. This updated formulation is known as the structured state space sequence layer (S4). Note that as a linear operator, the inverse discrete Fourier transform is amenable to backpropagation in the context of a neural network. Note also that under this formulation, the hidden state x(t) is only computed implicitly. Finally, [START_REF] Gu | On the parameterization and initialization of diagonal state space models[END_REF] presents a simplification of the S4 layer, known as diagonal S4 (S4D), in which A is approximated by a diagonal matrix.

Method

Given that the patch extraction of whole slide images at high magnifications results in long sequences of patches, we propose to incorporate a state space layer in a MIL aggregation network to better represent each patch sequence.

Neural state space models

In practice, neural state space models (SSM) simulate Equation 1 in discrete time, invoking a recurrence relation on the discretised hidden state,

x t = Ax t-1 + Bu t y t = Cx t + Du t (2)
where the sequences u t , x t , and y t are the discretised u(t), x(t), and y(t), and the modified model parameters arise from a bilinear discretisation [START_REF] Gu | Combining recurrent, convolutional, and continuous-time models with linear state space layers[END_REF]. As such, SSMs bear an inherent resemblance to RNNs, where the hidden representation x t can be interpreted as a memory cell for the observed sequence over the interval [0, t], and with Du t acting as a skip connection between the input and output at point t. Due to their lack of non-linearities, state space models can also be viewed as a convolution between two discrete sequences. Playing out the recurrence in Equation 2, one obtains,

y = K * u + Du, (3)
where u ∈ R L and y ∈ R L are the full input and output sequences, and the sequence K ∈ R L is defined as, which is computed efficiently by the S4D algorithm [START_REF] Gu | On the parameterization and initialization of diagonal state space models[END_REF]. Note that although SSM layers are linear, they may be combined with other, non-linear layers in a neural network. Note also that although Equation 3 is posed as modeling a onedimensional signal, in practice multi-dimensional inputs are modelled simply by stacking SSM layers together, followed by an affine "mixing" layer.

K = (CB, CAB, . . . , CA L-1 B), (4)

MIL training

In our pipeline (Figure 1) WSIs are first divided into a sequence of L patches {u 1 , u 2 , . . . , u L }, where L will vary by slide. A pretrained ResNet50 is then used to extract a 1024-dimensional feature vector from each patch {u 1 , u 2 , . . . , u L }, which constitute the model inputs. We define a SSM-based neural network F to predict a WSI-level class probability given this input sequence, ŷ = F({u 1 , u 2 , . . . , u L }).

(

) 5
The architecture of F is composed of an initial linear projection layer, used to lower the dimensionality of each vector in the input sequence. A SSM layer is then applied feature-wise by applying the S4D algorithm. That is, Equation 3, including the skip connection, transforms the sequence {u 1,d , u 2,d , . . . , u L,d } for all features d, and the resulting sequences are concatenated. A linear "mixing" layer is applied token-wise, doubling the dimensionality of each token, followed by a gated linear unit [START_REF] Dauphin | Language modeling with gated convolutional networks[END_REF] acting as an output gate, which restores the input dimensionality. For the SSM layer, we used the official implementation of S4D1 . A max pooling layer merges the SSM layer outputs into a single vector, which is projected by a final linear layer and softmax to give the class probabilities ŷ. The model is trained according to,

L M IL = - 1 M M m=1 log ŷcm , (6)
where ŷcm denotes the probability corresponding to c m , the slide-level label of the sequence corresponding to the m th of M whole slide images.

Multitask training

One advantage of processing an entire slide as a sequence is the ease with which additional supervision may be incorporated, when available. A patch-level ground truth creates the opportunity for multitask learning, which can enhance the representations learned for slide-level classification. As an extension of our base model in Equation 6, we train a multitask model to jointly predict a slidelevel and patch-level labels. Prior to the max pooling layer of the base model, an additional linear layer is applied to each sequence token, yielding L additional model outputs. This multitask model is trained according to a sum of log losses,

L M T = - 1 M M m=1 log ŷcm + λ L • L l=1 log ŷc m,l , (7)
where c m,l indexes the class of the l th patch in the m th training slide and λ is a tunable hyperparameter used to modulate the relative importance of each task.

Implementation details

We extracted patches of size 256 × 256 from the tissue regions of WSIs at 20x magnification. Following CLAM [START_REF] Lu | Data-efficient and weakly supervised computational pathology on whole-slide images[END_REF], the third residual block of a pretrained ResNet50 [START_REF] He | Deep residual learning for image recognition[END_REF] was used as a feature extractor, followed by a mean pooling operation, resulting in a 1024-dimensional representation for each patch. These features were used as inputs to all models. All model training was performed under a 10-fold cross-validation, and all reported results are averaged over the validation sets of the folds, aside from CAMELYON16, for which the predefined test set was utilized. Thus, for CAMELYON16, we report test set performances averaged over the validation.

Baseline models were chosen to be prior art CLAM [START_REF] Lu | Data-efficient and weakly supervised computational pathology on whole-slide images[END_REF] and TransMIL [START_REF] Shao | Transmil: Transformer based correlated multiple instance learning for whole slide image classification[END_REF]. The official code of these two models was used to perform the comparison. In addition, we included a vanilla transformer, a LSTM RNN, and based on mean and max pooling. Our vanilla transformer is composed of two stacked self-attention blocks, with four attention heads, a model dimension of 256, and a hidden dimension of 256. For the LSTM, we used an embedding size of 256 and a width of 256. The pooling models applied pooling feature-wise across each sequence, then used a random forest with 200 trees for classification. For the S4 models, the dimension of the state matrix A was tuned to 32 for CAMELYON16 and TCGA-RCC, and 128 for TCGA-LUAD. Our models were trained using the Adam [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] optimizer with the lookahead method [START_REF] Zhang | Lookahead optimizer: K steps forward, 1 step back[END_REF], with a learning rate of 2 • 10 -4 , and weight decay of 10 -4 for TCGA-LUAD and TCGA-RCC and 10 -3 for CAMELYON16. Early stopping with a patience of 10 was used for all our training. Our implementation is publicly available2 .

Experiments and discussion

4.1 Data CAMELYON16 [START_REF] Litjens | h&e-stained sentinel lymph node sections of breast cancer patients: the camelyon dataset[END_REF] is a dataset that consists of resections of lymph nodes, where each WSI is annotated with a binary label indicating the presence of tumour tissue in the slide, and all slides containing tumors have a pixel-level annotation indicating the metastatic region. In multitask experiments, we use this annotation to give each patch a label indicating local tumour presence. There are 270 WSIs in the training/validation set, and 130 WSIs in the predefined test set. In our experiments, the average patch sequence length arising from CAMELYON16 is 6129 (ranging from 127 to 27444).

TCGA-LUAD is a TCGA lung adenocarcinoma dataset that contains 541 WSIs along with genetic information about each patient. We obtained genetic information for this cohort using Xena browser [START_REF] Goldman | Visualizing and interpreting cancer genomics data via the xena platform[END_REF]. As a MIL task, we chose the task of predicting the patient mutation status of TP53, a tumor suppressor gene that is highly relevant in oncology studies. The average sequence length is 10557 (ranging from 85 to 34560).

TCGA-RCC is a TCGA dataset for three kidney cancer subtypes (denoted KICH, KIRC, and KIRP). It consists of 936 WSIs (121 KICH, 518 KIRC, and 297 KIRP). The average sequence length is 12234 (ranging from 319 to 62235).

Results

Multiple instance learning results. We evaluate our method on each dataset by accuracy and area under receiver operating characteristic curve (AUROC). For multiclass classification, these were computed in a one-versus-rest manner.

Table 1 summarises the comparison between our proposed model and baselines. For the CAMELYON16 dataset, our method performs on par with Trans-MIL and the CLAM models, while it clearly outperforms the other methods. Similarly, in the TCGA-LUAD dataset the proposed model achieves comparable performance with both CLAM models, while outperforming TransMIL and the other methods. We note that TCGA-LUAD proves to be a more challenging dataset for all models. Moreover, our method outperforms CLAM models on the TCGA-RCC dataset, while reporting very similar performance with respect to TransMIL. Overall, looking at the average metrics per model across all three datasets, our proposed method achieves the highest accuracy and the second highest AUROC, only behind CLAM-MB. A pairwise t-test between the proposed method, CLAM, and TransMIL shows that there is no statistical significance performance difference (see supplementary material).

We further compare our method with respect to model and time complexity. In Table 2 we report the number of trainable parameters, as well as the inference time for all models. The number of parameters is computed with all models configured to be binary classifiers, and the inference time is computed as the average time over 100 samples for processing a random sequence of 1024-dimensional vectors of length 30000. For our proposed method, we report both models with the [START_REF] Vaswani | Attention is all you need[END_REF] 1 054 978 2.60 LSTM [START_REF] Graves | Long short-term memory. Supervised sequence labelling with recurrent neural networks[END_REF] 789 250 320.52 CLAM SB [START_REF] Lu | Data-efficient and weakly supervised computational pathology on whole-slide images[END_REF] 790 791 0.84 CLAM MB [START_REF] Lu | Data-efficient and weakly supervised computational pathology on whole-slide images[END_REF] 791 048 5.85 TransMIL [START_REF] Shao | Transmil: Transformer based correlated multiple instance learning for whole slide image classification[END_REF] 2 Table 3 shows the effect of modifying parts of the architecture on the results for TCGA-RCC. Most modifications had very little impact on AUROC, but a more significant impact can be seen on the accuracy of the model. Models A and B show that stacking multiple SSM layers results in lower accuracy, which was observed over all three datasets, while models C and D show that modifying the state dimension of the SSM module can have an impact on the accuracy. The optimal state space dimension varies depending on the dataset.

Multitask learning results. We explored the ability of our model to combine slide-and patch-level information on the CAMEYLON16 dataset. We compared our model with the best performing model on CAMELYON16, TransMIL. Both models were trained according to Equation 7 with λ = 5 tuned by hand. In Table 4 we give slide-level accuracy and AUROC for the two models. We observe that all accuracies and AUROC increase compared with those reported in Table 1. This indicates that the use of patch-level annotations complements the learning of the slide-level label. We furthermore observe that our model outperforms TransMIL when combining slide-and patch-level annotations. We map the sequence of output probabilities to their slide coordinates giving a heatmap localising metastasis (see supplementary material).

Performance on longest sequences. In order to highlight the inherent ability of SSM models to effectively model long sequences, we performed an experiment on only the largest WSIs of the TCGA-RCC dataset. Indeed, this dataset contains particularly long sequences (up to 62235 patches at 20x). We evaluated the trained models for each fold on a subset of the validation set, only containing sequences with a length in the 85 th percentile. Table 5 shows the obtained average accuracy (weighted by the number of long sequences in each validation set) and AUROC on both CLAM models, TransMIL, and our proposed method. Both in terms of AUROC and accuracy, our method outperforms the other methods on long sequences, while the performances are comparable to Table 1, albeit slightly lower, illustrating the challenge of processing large WSIs.

Conclusions

In this work we have explored the ability of state space models to act as multiple instance learners on sequences of patches extracted from histopathology images. These models have been developed for their ability to memorise long sequences, and they have proven competitive with state of the art MIL models across a range of pathology problems. Additionally, we demonstrated the ability of these models to perform multiclass classification, which furthermore allowed us to visualise the localisation of metastasic regions. Finally, we demonstrated that on the longest sequences in our datasets, state space models offer better performance than competing models, confirming their power in modeling long-range dependencies.

Fig. 1 :

 1 Fig. 1: Overview of the proposed pipeline. In the first step, patches are extracted from a regular grid on a WSI. These patches are embedded using a pre-trained ResNet50 and are aggregated by a sequence model based on a state space layer.

Table 1 :

 1 Comparison of accuracy and AUROC on three datasets CAMELYON16, TCGA-LUAD, TCGA-RCC, and on average. All metrics in the table are the average of 10 runs. Best performing methods are indicated in bold and second best in italics. * indicates results from[START_REF] Shao | Transmil: Transformer based correlated multiple instance learning for whole slide image classification[END_REF].

	Dataset	CAMELYON16 TCGA-LUAD	TCGA-RCC	Average
	Metric	Acc. AUROC Acc. AUROC Acc. AUROC Acc. AUROC
	Mean-pooling 0.5969 0.5810 0.6261 0.6735 0.8608 0.9612 0.6946 0.7386
	Max-pooling 0.7078 0.7205 0.6328 0.6686 0.8803 0.9659 0.7403 0.7850
	Transformer [21] 0.5419 0.5202 0.5774 0.6214 0.7932 0.9147 0.6375 0.6854
	LSTM [8]	0.5310 0.5053 0.5389 0.5208 0.6654 0.7853 0.5784 0.6038
	CLAM SB [17] 0.8147 0.8382 0.6859 0.7459 0.8816 * 0.9723 * 0.7941 0.8532
	CLAM MB [17] 0.8264 0.8523 0.6901 0.7573 0.8966 * 0.9799 * 0.8044 0.8632
	TransMIL [19] 0.8287 0.8628 0.6348 0.7015 0.9466 * 0.9882 * 0.8034 0.8508
	Ours	0.8217 0.8485 0.6879 0.7304 0.9426 0.9885 0.8174 0.8558

Table 2 :

 2 Comparison of parameter count and inference time for all methods.

	Model	Number of parameters Inference time (ms)
	Mean-pooling	1 025	5.60
	Max-pooling	1 025	77.49
	Transformer

Table 3 :

 3 Ablation study for the different SSM components on the TCGA-RCC dataset. Best results in bold.

	Model SSM layers State dimension Accuracy AUROC
	A	2	32	0.9236 0.9813
	B	3	32	0.9179 0.9834
	C	1	16	0.9352 0.9846
	D	1	64	0.9352 0.9861
	Ours	1	32	0.9426 0.9885

Table 4 :

 4 Comparison of accuracy and AUROC for models trained as multitask classifiers on the CAMELYON16 dataset. Best results in bold.

	Model	Accuracy AUROC
	TransMIL [19] 0.8403 0.8828
	Ours	0.8488 0.8998

Table 5 :

 5 Results of CLAM SB, CLAM MB, TransMIL, and our proposed method on long sequences. Best results in bold.

	Model	Accuracy AUROC
	CLAM SB [17] 0.9149 0.9635
	CLAM MB [17] 0.8936 0.9654
	TransMIL [19] 0.9007 0.9652
	Ours	0.9220 0.9737

https://github.com/HazyResearch/state-spaces

https://github.com/MICS-Lab/s4 digital pathology

Acknowledgments This work has benefited from state financial aid, managed by the Agence Nationale de Recherche under the investment program integrated into France 2030, project reference ANR-21-RHUS-0003. This work was partially supported by the ANR Hagnodice ANR-21-CE45-0007. Experiments have been conducted using HPC resources from the "Mésocentre" computing center of CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France.


Leo Fillioux

Emilie Gontran

Jérôme Cartry

Jacques Rr Mathieu

Sabrina Bedja

Alice Boilève

Paul-Henry Cournède

Fanny Jaulin

Stergios Christodoulidis

Maria Vakalopoulou

Spatio-Temporal Analysis of Patient-Derived Organoid Videos Using Deep Learning for the Prediction of Drug Efficacy

Over the last ten years, Patient-Derived Organoids (PDOs) emerged as the most reliable technology to generate ex-vivo tumor avatars. PDOs retain the main characteristics of their original tumor, making them a system of choice for pre-clinical and clinical studies. In particular, PDOs are attracting interest in the field of Functional Precision Medicine (FPM), which is based upon an ex-vivo drug test in which living tumor cells (such as PDOs) from a specific patient are exposed to a panel of anti-cancer drugs. Currently, the Adenosine Triphosphate (ATP) based cell viability assay is the gold standard test to assess the sensitivity of PDOs to drugs. The readout is measured at the end of the assay from a global PDO population and therefore does not capture single PDO responses and does not provide time resolution of drug effect. To this end, in this study, we explore for the first time the use of powerful large foundation models for the automatic processing of PDO data. In particular, we propose a novel imaging-based high-throughput screening method to assess real-time drug efficacy from a time-lapse microscopy video of PDOs. The recently proposed SAM algorithm for segmentation and DI-NOv2 model are adapted in a comprehensive pipeline for processing PDO microscopy frames. Moreover, an attention mechanism is proposed for fusing temporal and spatial features in a multiple instance learning setting to predict ATP. We report better results than other non-time-resolved methods, indicating that the temporality of data is an important factor for the prediction of ATP. Extensive ablations shed light on optimizing the experimental setting and automating the prediction both in real-time and for forecasting.

Introduction

Precision medicine aims to optimize the choice of drug given the characteristics of the patient, so as to optimize certain aspects such as the efficacy of the treatment or qual-ity of life of the patient. Although doing this on a caseby-case basis by a clinician seems impractical, artificial intelligence-driven tools help guide this approach. In this objective, FPM [START_REF] Letai | Functional precision cancer medicine-moving beyond pure genomics[END_REF] bases this optimization on tests performed on live patient cells.

Patient-derived organoids (PDOs) have gained great interest over the last few years as they represent minimalistic models to mimic essential features from the tissue they originate from. In the context of cancer therapy, drug efficiency can be limited due to the development of resistance in patients as well as other evolutionary changes in the tumors over time. PDOs represent a good testbed for physicians, researchers, and patients to assess personalized drug efficacy, based on tumor-specific patient characteristics.

The gold standard method to assess drug efficacy on cells relies on Adenosine Triphosphate (ATP), an energy molecule released in active cells through metabolic reactions. ATP is a biomarker for evaluating the number of viable cells. The quantity of released ATP measured by luminescence is proportional to the number of living cells. Thus, ATP quantity assessed by luminescence counts serves as a readout to evaluate drug efficacy with the estimation of remaining living cells in the experimental sample. This test is easily implementable in bench labs and reliable for assessing the global cell population response to drug exposure [START_REF] Phan | A simple high-throughput approach identifies actionable drug sensitivities in patient-derived tumor organoids[END_REF]. However, since it causes cell lysis, the ATP test is a destructive assay, that does not allow to assess real-time organoid drug response, and post-ATP observations to evaluate long-term viability changes and drug resistance.

The emergence of large foundation models in various fields of machine learning has allowed for novel solutions in a variety of downstream tasks. Notably, SAM (Segment Anything Model) [22] or SEEM (Segment Everything Everywhere All at Once) [START_REF] Zou | Segment everything everywhere all at once[END_REF] have facilitated segmentation by providing a general model for segmentation from various types of prompts, which is extremely valuable for tasks for which little to no annotations can be found. Similarly, self-supervised [START_REF] Chen | A simple framework for contrastive learning of visual representations[END_REF][START_REF] Bardes | Vicreg: Variance-invariance-covariance regularization for selfsupervised learning[END_REF][START_REF] Goyal | Self-supervised pretraining of visual features in the wild[END_REF]9] or unsupervised models such as DINOv2 [START_REF] Oquab | Dinov2: Learning robust visual features without supervision[END_REF] provide high-quality descriptors for data and facilitate the training for downstream tasks with relatively small data and task-specific samples. This training setting allows the extracted features to be task-agnostic and, therefore, to adapt more easily to new tasks. Building on these recent advances, in this paper, we propose a new high-throughput screening method for the analysis of PDOs testing multiple drugs and high-quality videos. Indeed, processing of this data is usually performed in a manual and time-consuming setting.

In this paper, our contributions are the following. To our knowledge, we present the first fully automatic method for processing high-quality videos of PDOs, conducting a spatio-temporal analysis for the prediction of ATP. We propose an efficient, automatic, and accurate prompt engineering paradigm taking into account the temporal characteristics of PDOs for using SAM without the need for additional training. Finally, we explore powerful recent foundation models for the spatio-temporal representation of PDOs for the first time coupling them with time sequence modeling and multiple instance learning.

An extensive experimental analysis was performed to identify the best representations for PDOs as well as the most informative time frames, which we used for the accurate prediction of ATP. This opens the possibility to integrate any other clinical endpoint, to match PDO drug sensitivity and patient clinical response.

Related work

Foundation Models in Computer Vision. Foundation models are precious resources and stimulate research by both providing a strong solution for a specific task and by proposing a novel innovative approach. The authors of the "Florence" model [START_REF] Yuan | Florence: A new foundation model for computer vision[END_REF] introduced a method for learning joint visual-textual features that can be adapted to a multitude of joint tasks. The CLIP model [START_REF] Radford | Learning transferable visual models from natural language supervision[END_REF] uses a contrastive learning approach in order to learn image embeddings, guided by an associated textual description of the image, which can transfer zero-shot to a wide variety of tasks. This approach has led to many extensions [START_REF] Harold | Grounded language-image pre-training[END_REF][START_REF] Mu | Slip: Self-supervision meets language-image pretraining[END_REF]24]. DINOv2 [START_REF] Oquab | Dinov2: Learning robust visual features without supervision[END_REF] is a deep learning model trained in a selfsupervised manner that leverages more stable training at a bigger scale for better taking advantage of large datasets. The model is trained on a dataset of 142M images. It is based on the ViT architecture [START_REF] Dosovitskiy | An image is worth 16x16 words: Transformers for image recognition at scale[END_REF] and is made available in different sizes. Very recently, foundation models for image segmentation such as SegGPT [START_REF] Wang | Seggpt: Segmenting everything in context[END_REF], SEEM [START_REF] Zou | Segment everything everywhere all at once[END_REF], and SAM [22] have surfaced. SAM provides a strong zeroshot transfer approach by combining a powerful image encoder with a prompt encoder which can adapt to multimodal prompts. Its convincing adaptability led to many adaptations for videos [START_REF] Yang | Track anything: Segment anything meets videos[END_REF], lightweight versions [49], or medical applications [START_REF] Wu | Medical sam adapter: Adapting segment anything model for medical image segmentation[END_REF][START_REF] Liu | Samm (segment any medical model): A 3d slicer integration to sam[END_REF]. While these models show impressive results and generalizability, their adaptations to medical ap-plications often show limitations. For example, recent generative models such as GLIDE [START_REF] Nichol | Glide: Towards photorealistic image generation and editing with text-guided diffusion models[END_REF] or DALL-E [START_REF] Ramesh | Zero-shot text-to-image generation[END_REF] show impressive results on natural images but struggle to generate realistic medical images without prior fine-tuning [21,2].

Time Sequence Modeling. Modeling the temporal aspect of data can be achieved either on extracted features of the frames or directly in an end-to-end fashion on the videos. Once features have been extracted from frames of the videos, standard sequence models such as LSTMs [START_REF] Hochreiter | Long short-term memory[END_REF] or Transformers [START_REF] Vaswani | Attention is all you need[END_REF][START_REF] Zerveas | A transformerbased framework for multivariate time series representation learning[END_REF] can be applied for a variety of tasks. They aim to capture dependencies between elements of a sequence, either in the form of recurrent neural networks or self-attention mechanisms. State space modeling is gaining attraction in sequence modeling, especially for long sequences [START_REF] Gu | Efficiently modeling long sequences with structured state spaces[END_REF], and more specifically for modeling time series [START_REF] Zhang | Effectively modeling time series with simple discrete state spaces[END_REF]. For end-to-end video modeling, Transformer adaptations are very popular such as the Video Swin Transformer [27], ViViT [START_REF] Arnab | Vivit: A video vision transformer[END_REF], or the TimeSformer [START_REF] Bertasius | Is space-time attention all you need for video understanding?[END_REF]. We are also beginning to see adaptations of large language models (LLMs) for video understanding tasks [START_REF] Chen | Videollm: Modeling video sequence with large language models[END_REF]. However, even if these models are available, their application to PDOs has not yet been explored and investigated, especially in a low-availability data regime.

Analysis of PDO Microscopy Images. With the emergence of PDOs as a key technology for FPM in the last few years, machine learning methods to analyze this data have followed. Earlier studies [START_REF] Bian | A deep learning model for detection and tracking in high-throughput images of organoid[END_REF][START_REF] Matthews | Organoid: a versatile deep learning platform for organoid image analysis[END_REF] looked for methods for tracking and analyzing the dynamics of PDOs through time. Authors of "D-CryptO" [START_REF] Abdul | D-crypto: deep learning-based analysis of colon organoid morphology from brightfield images[END_REF] proposed to classify types of organoids based on their morphology. Later research sought to use organoid analysis pipelines for other downstream tasks, such as the prediction of kidney differentiation [START_REF] Park | Deep learning predicts the differentiation of kidney organoids derived from human induced pluripotent stem cells[END_REF] or the prediction of a biomarker of Huntington's disease [START_REF] Metzger | Deep-learning analysis of micropatternbased organoids enables high-throughput drug screening of huntington's disease models[END_REF]. Recently [7] have developed a pipeline for predicting ATP from microscopy images at single time points in a multiple instance learning setting. However, all these methods focus on analyzing single representations of organoids, losing a lot of information that arises from their temporal dynamics for different drug treatments. Moreover, they are mainly based on classical processing techniques, such as the extraction of predefined visual features, the extraction of ResNet features, or segmentation through biological staining.

Methods

In this section, we will start by describing the specificities of the dataset and the preprocessing steps, and then proceed to an in-depth presentation of the proposed method. Specifically, we will present the methodology for the segmentation of the cavities and organoids by employing SAM [22], then we will present the details of the feature extraction from the organoid regions using DINOv2 [START_REF] Oquab | Dinov2: Learning robust visual features without supervision[END_REF]. The model for the final prediction of the ATP given the features from the cavities under the multiple instance learning framework will then be presented. Figure 1 illustrates the proposed pipeline.

Data

To capture real-time information about drug efficacy on patient organoids, bright-field imaging techniques represent a good non-invasive alternative to the ATP bioluminescence test. However, having single-organoid resolution may be a challenging task when working with the standard experimental setup where organoids are embedded in hydrogel droplets, which results in heterogenous organoid distribution at different depths and subsequently heterogeneous organoid sizes, making the drug efficacy assessment less accurate. To overcome these issues, we developed a highthroughput experimental and imaging pipeline. It consists in using high-throughput cell culture systems with 96-well plates containing 500 µm diameter cavities. Each cavity contains one organoid, while the setting allows for the parallel follow-up of individual organoids with a homogeneous depth distribution. Testing multiple drugs in parallel thus becomes easier with one image containing the information of multiple organoids under the same treatment. One type of drug at a given concentration is applied to each well, and the ATP can only be measured on the well level. Figure 2 shows one well containing multiple cavities together with their automatic segmentation. A total of 116 wells are imaged, from which we detect 8241 cavities, and therefore 8241 organoids.

A cavity is imaged in a set of frames T = {T 1 , ..., T f }, T i ∈ R w×h , where w and h are the frame dimensions, and f the number of frames in each timelapse. In this study, we developed a novel pipeline to segment each cavity S = {S 1 , ..., S f }, S i ∈ R w×h and using S and T , we extract interesting features from each organoid. Finally, we model a well by a set of cavity features C = {c 1 , . . . , c n } ∈ R n×k×f , where n is the number of cavities in the well, and k the dimensionality of the feature space. The final ATP value is defined as y ∈ R.

The patient-derived organoids (from now on called organoids, for simplicity) are extracted from colorectal cancer patients. In this analysis, we do not consider the concentration of the drugs and only observe their effects indirectly with the ATP. Drugs are introduced one day after the organoids have been formed in the cavities. The wells are imaged every 30 minutes for 100 hours, resulting in f = 200 frames for each well. Wells are digitized using an Agilent Lionheart FX digital microscope. Note that some experiments could take slightly longer than 200 frames, for those cases, only the last 200 frames of the experiments were retained, in order to keep the measurement of the ATP at the last frame.

Preprocessing.

Each scan consists of acquiring quarters of wells at three different depth levels, resulting in 12 images for each time point. A z-projection based on Sobel filters [START_REF] Kanopoulos | Design of an image edge detection filter using the sobel operator[END_REF] is used to project the images into the same plane, before stitching all four corners into a single image. Artifacts appear through time (e.g. evaporation of the liquid impacting the contrast), which may impact the segmentation pipeline. To account for this, we normalize the contrast of the timelapse through time. Finally, all frames of a timelapse are then coregistered using SIFT-based features [START_REF] David | Distinctive image features from scaleinvariant keypoints[END_REF]. All preprocessing is performed using ImageJ [START_REF] Schneider | Nih image to imagej: 25 years of image analysis[END_REF]. Figures 2a-2d show all steps of the preprocessing from the 12 input images per frame to the normalized frame.

Automatic Segmentation of Regions of Interest

In this study, we propose and develop an automatic pipeline based on SAM [22] for the automatic segmentation of regions of interest in the video of the organoids. For accurate segmentation, various characteristics of these regions are used and different prompt engineering strategies have been developed to handle the nature of the data and their temporal information.

Cavity Segmentation. The objective is to segment the individual cavities contained in the well, in order to represent the well as a set of cavities. We assume that the detection of the cavities only needs to be run in a single frame, as the frames of the video have been co-registered. The well is considered as the mask with the biggest surface area containing the center pixel in the frame is considered, as it is the biggest object in the frame. Given the candidate masks M = {m 1 , . . . , m n }, m well = m argmax i {area(mi)} . The cavities are defined as other regions detected by SAM, which have a circularity above a certain threshold and which are included in the detected well. For each point b i on the boundary B of a detected region m, we define P = {max i (dist(b i , B))}, which represents for each point on the contour, the distance to the point on the contour that is furthest away. The ratio min P max P defines the circularity. On a perfect circle, all values in P should be equal to the circle's diameter and therefore have circularity = 1. Once the cavities have been detected on a single frame, a crop around the cavity is extracted for all frames, giving a few dozen cavity timelapses per well. Figure 2e shows an example of the cavity segmentation masks on a normalized frame.

Organoid Segmentation. Segmentation of organoids is performed on the extracted videos of the cavities by designing the proper prompts for SAM [22]. The first prompt is generated for the last frame of the timelapse T f (as it is the frame where the organoid is the largest, and it is, therefore, more likely to find a point inside the region of interest). Canny edge detection [START_REF] Canny | A computational approach to edge detection[END_REF] is used to extract rough contours from the last frame, after having subtracted the average frame over time, in order to remove the background. These rough contours are then filtered using a series of binary morphological operators, before extracting a centroid, which is used as a positive point prompt for SAM. Figure 3 shows an example of these steps for the last frame. Once the mask for the last frame has been generated, the prompts will be generated backward from the last frame to the first frame. The center of the organoid mask is extracted from the prediction and an exponentially weighted average from the prompts from frames T t+1 to T t+10 is used to generate the prompt for frame T t . Multimask output is used to generate the top 3 masks, and the mask which has the highest Dice score with the generated mask from the previous frame (T t+1 ) is chosen. Postprocessing is performed on the mask by only selecting the detected region with the highest surface area. Algorithm 1 represents the whole organoid segmentation pipeline.

Inputs:

frames T = {T 1 , ..., T f }, T i ∈ R w×h Outputs:

segmentations S = {S 1 , ..., S f }, S i ∈ R w×h estimated boundary= Canny(T fmean(T )) prompt f = centroid(estimated boundary) S f = SAM(prompt f, T f ) all prompts = {prompt f} for i=(f-1); i=1; i --do prompt i = exp weighted average(all prompts, 10) S i = SAM(prompt i, T i ) all prompts.insert(prompt i) end for

Feature Extraction

Features are extracted from each frame of the cavity using the DINOv2 [START_REF] Oquab | Dinov2: Learning robust visual features without supervision[END_REF] model, which produces task-agnostic visual features from images. A lighter version (86M parameters) of the model is used, which produces k-dimensional feature vectors. Features are extracted per frame, using a crop of the cavity around the organoid, and by masking everything outside of the regions of interest, giving c i = DINOv2(T i ⊙ S i ). The model generates feature vectors of dimension k = 768 for each time frame.

ATP Prediction

Having access to features for each individual cavity but only having the ATP measure on the well level, the multiple instance learning setting seems to be the best approach to the problem. Each set of cavities C is associated with an ATP measure y. The ATP prediction model M should map an input C ∈ R n×k×f to an output ŷ ∈ R. The aggregation of the cavity representations on the well level is done by mean pooling in order to obtain a well-level representation mean(C) ∈ R k×f .

Each time frame may have a different impact on the prediction of the ATP. In fact, we expect later time frames to have more impact than earlier time frames. We give the model the freedom to learn how much weight to put to each time frame using a normalized weight vector w t ∈ R f that is used in a weighted average across time frames. This ensures that all cavities across all wells share the same importance for a given time frame. A feature-wise weight vector w k ∈ R k is used in order to learn the relative importance of each feature that is shared across all wells.

The last element of the model is a multilayer perceptron (MLP) which is simply composed of four linear layers, with PReLU [START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[END_REF] activation function.

M : R n×k×f → R C → MLP(mean(C) • w t ∥w t ∥ ⊙ w k ) (1)
The loss function used for training is composed of two elements, a relative L 1 loss, to directly optimize for the mean absolute percentage error (MAPE), and a L 1 loss, which has been normalized with the maximum ATP value in the training set to have a comparable range to the relative L 1 loss.

L rel (y, ŷ) = |y -ŷ| y (2) L norm (y, ŷ, y max ) = |y -ŷ| y max (3)
L(y, ŷ, y max ) = αL rel + (1 -α)L norm ( 4
)
Where y is the ATP ground truth, ŷ the output of our model, y max the maximum value of ATP in the training set, and α the weight given to L rel . A relative L 1 loss is used because of the varying range that the ATP can take within experiments. The normalized L 1 loss is used to ensure that the model still has an acceptable performance on higher values of ATPs.

Implementation Details

Over all experiments, training was performed under the same scheme. Cross-validation is performed over 4 splits, on the well level, resulting in training sets of size 87 and validation sets of size 29. The reported results are the average over all folds of the validation sets. All trainings were performed using the PyTorch deep learning library [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] in Python. The AdamW [START_REF] Loshchilov | Decoupled weight decay regularization[END_REF] optimizer was used with a learning rate of 1 • 10 -3 , weight decay of 1 • 10 -1 , and a batch size of 8. The model was trained over 2000 epochs, and early stopping with a patience of 200 epochs was used. The α parameter which balances both terms of the loss was set at α = 0.5 after grid search. All trainings were accelerated using Nvidia Tesla V100 GPUs.

Experiments

To evaluate the performance of the method we used the MAPE and Pearson correlation coefficient metrics. All performances are measured on each validation set, utilizing all four models trained on the corresponding training sets. The MAPE, measures the mean relative error between the prediction ŷ and the label y. Because of the very wide range of values of the ATP (both within an experiment and between different experiments) as presented in Table 1, the MAPE appears suitable for ensuring comparability, in particular for ATP bin ( 10

MAPE(y, ŷ) =

|y -ŷ| y

The Pearson correlation coefficient quantifies the linear correlation between the set of predictions Ŷ and the set of labels Y . Higher values of the Pearson correlation coefficient indicate better performance.

Pearson(Y, Ŷ ) = i (Y i -Ȳ )( Ŷi -Ȳ ) i (Y i -Ȳ ) 2 ( Ŷi -Ȳ ) 2 (6)
Figure 4 qualitatively compares three methods of organoid segmentation to justify the superiority of our proposed method over simple intensity-based methods (Otsu thresholding) and justifies the use of our prompt generation for the use of SAM [22].

Results

With our proposed high-throughput screening method, we obtain an average MAPE and Pearson correlation coefficient on the four validation sets of 0.1755 and 0.9214, respectively. The relative error of 17.55% is to be put in context with the wide range of ATP values in our dataset, with a ratio between the highest and lowest values of ATP of 14. Table 1 shows the distribution of ATP values along with the corresponding performance of the model in terms of MAPE in each bin. The model performs worst on the bin with the smallest values of ATP, which may be explained by the precision of the measurement of ATP.

Feature Extraction. In Table 2, we explore the impact of four different feature extraction methods on the final results. Different classical and deep learning features have been explored to highlight the best representations for organoids. The pyradiomics Python package [START_REF] Joost | Computational radiomics system to decode the radiographic phenotype[END_REF] allows the extraction of imaging features from regions of interest. For this study, we extracted features related to first-order statistics, gray-level matrices, and shape features. A total of 93 features have been extracted per frame. ResNet50 [START_REF] He | Deep residual learning for image recognition[END_REF] (pre-trained on ImageNet [12]) is a popular classification architecture, which can be used for feature extraction when removing the classification head. It produces 2048dimensional features. VICReg [START_REF] Bardes | Vicreg: Variance-invariance-covariance regularization for selfsupervised learning[END_REF] is a feature extraction model trained in a self-supervised manner. For our study, we used a model based on the ResNet50 architecture, which produces 4096-dimensional features. Features extracted using models trained in a self-supervised manner (VICReg and DINOv2) seem to adapt better to different downstream tasks, compared to models trained for other specific tasks (ResNet50 on ImageNet), or classical predefined features. The superiority of our proposed model is highlighted both in terms of MAPE and correlation coefficient.

Attention Mechanism. We used two types of attention mechanisms, one attending to the temporal aspect of the data, and the other attending to the individual features. In Table 3 we perform an ablation study for three variations for each type of attention mechanism. For the temporal attention, we explored the impact of only using the last frame, multihead attention [START_REF] Vaswani | Attention is all you need[END_REF] (MHA), or the learnable weight vector w t . Similarly, for the feature-wise attention, we explored the impact of having no such attention, using multihead attention (MHA) or the learnable weight vector w k . In terms of temporal attention, using w t clearly outperforms other methods, especially the multihead attention, while using the last frame still gives reasonable results. w t seems to grasp the temporal importance of each frame. Concerning the feature-wise attention, the use of multihead attention clearly does not seem adequate for this application, while the effect of incorporating w k compared to its absence is not readily apparent. The analysis of its weights could however be used for finer explainability of the model. Cavity Aggregation. In multiple instance learning, the choice of the bag level aggregation plays a crucial role. In Table 4 we tested multiple methods of aggregation, minpooling, max-pooling, mean, sum (which is different from the mean as the number of cavities per bag varies), and using a SE [19] block followed by a sum aggregator. The sum and mean operators outperform other types of aggregation, which can be explained by the fact that the ATP is a direct function of the number of organoid cells in the well. This information is not captured by the min-pooling and maxpooling operators. The SE block learns the weight to give to each cavity based on its features, which is used to model the fact that different cavities contribute differently to the ATP. While this may be true, we hypothesize that the information is already present in the features, and the use of the SE block therefore only adds noise.

Takeaway. Our model has a mean MAPE 0.1755 and Pearson of 0.9214 on the validation set. Trying out different feature extraction methods shows that features from DI-NOv2 [START_REF] Oquab | Dinov2: Learning robust visual features without supervision[END_REF] provide the best results. Comparing the impact of a variety of feature and temporal attention schemes justifies the use of w t and w k , while we also show the superiority of using the mean as the bag-level aggregation function.

Learned Attention Weights

Once the model is trained, analyzing the weights learned by the w t and w k vectors gives insights into our experiments. Figure 5a shows the weight associated with each frame in the w t vector after training. Intuitively, we expect later frames to have the most importance. Without any constraints during training, we see that the learned temporal attention has learned a continuous and smooth distribution, meaning that nearby frames will have relatively similar importance. It is noteworthy that the initial frames seem to have more importance for the final prediction than frames after 24 hours of the experiment. This could be attributed to the fact that the organoid's initial state is significantly determinant of the final ATP, likely due to its size, as larger organoids tend to exhibit higher values of ATP.

Examining the feature importance in Figure 5b shows that the feature attention does not highlight any specific feature, as indicated by a ratio of 2 between the feature with the highest and lowest attention weight. Note that in Figure 5b, the features are sorted by their attention weight.

Takeaway. The temporal attention vector w t learns meaningful relative importance of the time frames, while the feature-wise attention vector w k does not make any particular feature stand out.

Comparison to SOTA

To the best of our knowledge, only Ins-ATP [7] have explored the use of machine learning for the prediction of ATP from organoid microscopy images. They imaged multiple organoids placed in a single matrigel drop (compared to our method which has one organoid per cavity) at a single time point, which yields very different visual results to our images. We implemented the Ins-ATP method and adapted it to fit our data, by defining the bags of instances as a set of cavities from the wells, only using the last frame. This can be thought of as a compromise between the presented MeshIns-ATP and DeepIns-ATP [7] as the instances composing the bags are not learned by the model, but are chosen as meaningful regions of interest (the cavities). With this method, we obtained an average MAPE and Pearson correlation coefficient of the 4 validation sets of 0.4375 and 0.6794 respectively (compared to 0.1755 and 0.9214 respectively for our pipeline). This highlights the advances of our study in the development of a high-throughput screening method to assess real-time drug efficacy from a timelapse microscopy video of PDOs.

Takeaway. The Ins-ATP [7] method applied to our dataset gives lower performance than our proposed method.

Beyond Predicting Current ATP.

We have shown the ability of our proposed pipeline to, given a sequence of frames, predict the ATP measured at the last frame. Two main questions arise. How well does our model perform in forecasting the ATP? How many frames of history does our model need to predict the current ATP? Note that while we are retraining the model with features from a subset of frames, the segmentation maps are still the ones computed on all f frames.

Forecasting the ATP. Given our current data (i.e. a sequence of f = 200 frames, and a measure of ATP at the last frame), we can evaluate the performance of our pipeline for forecasting the ATP by training the model while omitting the last frames. In Figure 6a, we show the performance of the model when trained only on frames from 0 to i, equivalent to predicting the ATP (f -i) frames in advance. As expected, the performance of the model drops as the model is trained to predict ATPs earlier in the video. However, predicting 2 days in advance a MAPE of 0.28 seems reasonable compared to the best performance of 0.1755, especially when considering the wide range of ATP values.

Required history for ATP prediction. Similarly, if we omit the first frames from the training, we can evaluate how many frames of history our model needs for predicting the ATP. In Figure 6b we show the performance of the model when trained only on frames from i to f (i.e. having access to f -i frames). The drop in performance is not as clear as the one shown in Figure 6a. This indicates that our model seems to perform well when not taking too many frames into account: the peak performance appears to happen when about 15 frames are given to the model. This is positive because it suggests that our model could be used for the live measurement of ATP as the images are still being acquired, as it does not need a full history of the organoid cells. However, it is important to note that although the ATP prediction part of the model seems to perform well with only 15 frames, the segmentation pipeline (and more precisely, the generation of the first prompt for SAM) relies on movement within the timelapse. We tested how well the segmentation of the organoids worked on the last 15 frames by computing the Dice score with the last 15 frames of the segmentation map computed using all frames as ground truth. The mean Dice score is 0.80, with a very uneven distribution among cavities: 75% of the cavities obtain a Dice over 0.90, while 15% obtain a Dice under 0.20. It seems that most segmentations are not affected by the use of a smaller number of frames, but for those cases which are affected, the segmentations are prone to complete failure.

Takeaway. Our model can predict the ATP in advance with a reasonable drop in performance and is able to predict with only about 15 frames of history making it usable for online predictions.

Conclusion

The estimation of the ATP is a standard method for the estimation of drug efficacy on organoids. However, it allows the measurement at a single timepoint for the entire experiment. In this paper, we propose a method for the spatiotemporal analysis of organoid microscopy timelapse videos, for the prediction of ATP. We assess the performance of our approach for predicting the current ATP with different ablations and we report better performance than SOTA.

Future work includes the further exploration of foundation models for the analysis of organoid videos. In this study, the foundation models are used as frozen blocks. However, authors of DINOv2 [START_REF] Oquab | Dinov2: Learning robust visual features without supervision[END_REF] showed that finetuning the model encoders to a specific dataset improved the results on the dataset. One step could be to finetune the feature extraction model to improve the representations of the organoids to be better adapted to the prediction of ATP. Similarly for the SAM [22] model, we could finetune the mask decoder part of the model, or learn more efficient prompts.

Figure 1 .

 1 Figure 1. Illustration of our proposed pipeline. (a) The overall pipeline taking as input the whole well timelapse and outputting the predicted ATP ŷ. (b) Organoid segmentation based on automatic generation of prompts for SAM [22]. (c) Feature extraction given input frames and associated segmentation maps focusing in the region of interest and using DINOv2 model [34].

  Figure 2. Example of the steps needed for preprocessing and cavity segmentation. (a) All images needed to construct one frame (4 corners at 3 different z-levels each). (b) Projected corners. (c) Stitched frame. (d) Local contrast normalization applied to the stitched frame. (e) Cavity segmentation masks.

Figure 3 .Algorithm 1

 31 Example of the steps used in the generation of the prompt for the segmentation of the last frame of a cavity timelapse. (a) Last frame of the timelapse. (b) Mean frame across time. (c) Difference between the last frame and the mean frame. (d) Result of Canny edge detection on (c). (e) Generated prompt for the frame. (f) Predicted mask. Organoid segmentation.

Figure 4 .

 4 Comparison of three methods of segmentation. (a) Our proposed method. (b) Otsu thresholding based segmentation. (c) Using SAM [22] without any prompts and choosing the mask with the highest predicted IOU. future studies which may use datasets with different distributions of ATP. Lower values of MAPE indicate better performance.

Figure 5 .

 5 Figure 5. Visualization of the temporal attention (a) and (sorted) feature attention (b) weights for a model after training (visualization for one fold only). Polynomial fit in red.

Figure 6 .

 6 MAPE (left) and Pearson correlation coefficient (right) as a function of the number of frames given to the model. (a) Performance of models trained only on frames from 0 to i. (b) Performance of models trained only on frames from i to 200.

Table 1 .

 1 5 ) [1.1, 2.8) [2.8, 4.5) [4.5, 6.2) [6.2, 8.0) [8.0, 9.7) [9.7, 11.4) [11.4, 13.1) [13.1, 14.9] Distribution of ATP values in our dataset, with associated MAPE.

	Count	23	19	11	15	0	7	17	24
	MAPE ↓	0.34	0.15	0.08	0.07	NA	0.10	0.23	0.12

Table 2 .

 2 Comparison of MAPE and Pearson correlation coefficient for four methods of feature extraction. Best results are indicated in bold, and second best results are underlined.

	Features Classical ResNet50 [16] VICReg [4] DINOv2 [34]
	MAPE ↓	0.3173	0.2174	0.1847	0.1755
	Pearson ↑	0.7883	0.8960	0.8939	0.9214

Table 3 .

 3 Comparison of MAPE and Pearson correlation coefficient for three types of temporal attention and feature-wise attention. Best results are indicated in bold, and second best results are underlined.

	Cavity aggregation	Min	SE [19]	Max	Sum	Mean
	MAPE ↓	0.2425 0.2138 0.2092 0.1797 0.1755
	Pearson ↑	0.8912 0.8922 0.9185 0.9204 0.9214

Table 4 .

 4 Comparison of MAPE and Pearson correlation coefficient for four methods of feature extraction. Best results are indicated in bold, and second best results are underlined.

Acknowledgments. This work has benefited from state financial aid, managed by the Agence Nationale de Recherche under the investment program integrated into France 2030, project reference ANR-21-RHUS-0003. Experiments have been conducted using HPC resources from the "Mésocentre" computing center.