Predicting Age-related Macular Degeneration Progression from Retinal Optical Coherence Tomography with Intra-Subject Temporal Consistency

Authors: arunava  chakravarty
taha  emre
dmitrii  lachinov
antoine  rivail
ursula  schmidt-erfurth
hrvoje  bogunović
Venue: MIDL_io_2024_Conference
Published: 2024-01-31 10:30:37

Abstract: The wide variability in the progression rates of Age-Related Macular Degeneration (AMD) and the absence of well-established clinical biomarkers make it difficult to predict an individual's risk of AMD progression from intermediate stage (iAMD) to late dry stage (dAMD) using Optical Coherence Tomography (OCT) scans.
To address this challenge, we propose to jointly train an AMD stage classifier to discriminate between iAMD and dAMD with a Neural-ODE that models the future trajectory of the disease progression in the learned embedding space. A temporal ordering is imposed such that the distance of a scan from the decision hyperplane of the AMD stage classifier is inversely related to its time-to-conversion. In addition, an intra-subject temporal consistency in the predicted conversion risk scores is ensured by incorporating a pair of longitudinal scans from the same eye during training. We evaluated our proposed method on a longitudinal dataset comprising 235 eyes (3,534 OCT scans) with 40 converters. The results demonstrate the effectiveness of our approach, achieving an average area under the ROC of 0.84 for predicting conversion within the next 6, 12, 18 and 24 months. Additionally, the Concordance Index of 0.78 surpasses the performance of several popular methods for survival analysis.

Main paper: 
 introduction: 
 Age-related macular degeneration (AMD) is a leading cause of blindness among the elderly population (Wong et al., 2014) . It is asymptomatic in its early and intermediate stages (iAMD), characterized by the presence of drusen. AMD gradually advances to the late stage leading to irreversible vision loss which could be categorized as either neovascular (nAMD) or dry (dAMD) . nAMD is caused by abnormal blood vessel growth in the choroid that leaks fluid into the retina. dAMD is more prevalent than nAMD and characterized by Geographic Atrophy (GA) due to the loss of Retinal Pigment Epithelium (RPE). Recently, for the first time, drugs for dAMD (Khanani et al., 2023; Heier et al., 2023) were approved by FDA. Patients in the iAMD stage are regularly monitored with longitudinal Optical Coherence Tomography (OCT) imaging across multiple visits to initiate treatment at the earliest onset of late AMD to minimize vision loss. Identifying iAMD patients at a high risk of dAMD conversion enables ophthalmologists to prioritize these cases for enhanced monitoring, facilitating early detection of dAMD onset. However, this is a challenging task due to the absence of well-established clinical biomarkers and significant inter-subject variations in the rate of AMD progression. Deep learning (DL) methods to predict the future risk of conversion of an eye from iAMD to dAMD can play a critical clinical role in supporting personalized treatments and clinical research by categorizing iAMD patients into distinct risk levels for biomarker identification and recruitment in clinical trials.Related Work: Existing methods for predicting the risk of conversion from iAMD to nAMD or dAMD fall into two main categories: biomarker and image-based approaches. Biomarker-based methods (Sleiman et al., 2017; Schmidt-Erfurth et al., 2018; Banerjee et al., 2020; de Sisternes et al., 2014; Lad et al., 2022) involve segmenting retinal tissues and pathologies to extract features, subsequently combined with clinical and demographic data for risk prediction. Notably, Banerjee et al. (2020) utilizes multiple past visit biomarkers in an LSTM network for future risk assessment. Image-based methods, however, directly utilize DL models on raw OCT scans, bypassing manual segmentation. A hybrid approach using both biomarker and image features for predicting nAMD conversion is presented in Yim et al. (2020) , employing an ensemble DL model. Unlabeled longitudinal OCT datasets have been used in Emre et al. (2022) ; Rivail et al. (2019) for feature learning via temporal self-supervised learning. These methods typically employ a binary classifier for predicting conversion within specific timeframes, (e.g., 2 years (Russakoff et al., 2019 ), 6 months (Yim et al., 2020; Emre et al., 2022) ), or multi-label classification for various discrete timeintervals (e.g., 6, 12, and 18 months (Rivail et al., 2019) ).The binary classification based approaches are limited by discretization of the conversion time and their inability to manage censoring, which occurs when an eye's actual conversion time is unknown due to missing follow-ups or non-conversion within a limited study duration. Survival analysis addresses these challenges. Discrete survival models are similar to multi-label classification but modify training loss to incorporate censoring and have recently been applied to predict dAMD conversion (Rivail et al., 2023) . A transformer model has also been used for discrete-time modeling of the hazard function from tabular clinical and demographic data (Hu et al., 2021) . Traditional non-DL continuous models of survival analysis have also been explored to capture AMD progression with handcrafted biomarkers using the linear Cox Proportional Hazard model (CoxPH) (Schmidt-Erfurth et al., 2018) . Although CoxPH has been extended with DL using images (Katzman et al., 2018), they have not yet been explored to model AMD progression so far. Moreover, these models are inflexible as each patient's hazard function is constrained to be a scaled version of the same baseline hazard across the entire population. SODEN (Tang et al., 2022) overcomes this issue by employing a N-ODE to model the cumulative hazard function for survival on tabular data. The GRU-ODE-Bayes (De Brouwer et al., 2019) proposed a N-ODE to extend the GRU based Recurrent Neural Network in continuous time, used in predicting disability progression in Multiple Sclerosis patients from tabular data of past history. Recently, N-ODEs have also been used to model the spatial evolution of GA segmentation in OCT (Lachinov et al., 2023) and Diabetic retinopathy in fundus images (Zeghlache et al., 2023) .Our key contributions are: (i) The time-to-conversion from iAMD to nAMD is modeled in continuous time, rather than discrete time-intervals as used in most existing methods. Our model can therefore use actual continuous conversion times as ground-truths during training and also predict conversion probabilities within arbitrary continuous times. (ii) Our novel N-ODE based modeling directly models the Cumulative Distribution Function(CDF) of the future conversion time instead of the cumulative hazard function used in existing methods like SODEN. Our SMGRU-ODE architecture also extends ODE-GRU by stacking multiple layers with multiple parallel heads. (iii) We incorporate intra-subject consistency by requiring the N-ODE estimates of the feature and risk at future time-points to be consistent with the values obtained using the actual OCT scan of the future visit. (iv) We jointly train a linear AMD stage classifier and employ a rank loss on its logits which is sensitive to censoring, to regularize the feature embedding. This facilitates patient stratification into risk groups based on a scalar risk score derived from the decision hyperplane distance for clinical studies or personalized treatment. Given an input OCT scan of an iAMD patient, the proposed method projects it to a feature embedding where the current feature is evolved over time with a N-ODE to forecast the future trajectory of the disease progression. The features estimated for any (continuous) future time-point are fed to a linear AMD stage classifier to predict the probability of the eye to have already converted within that time, thereby modeling the CDF of the future conversion time. In Survival analysis, the Ground Truth (GT) label for an OCT image I j is defined by the tuple (E j , T j ). The event indicator E j = 1 signifies that the eye associated with scan I j will progress from iAMD to dAMD, while E j = 0 denotes no conversion within the monitoring period. T j denotes the time of conversion from the current visit if E j = 1 or the censoring time until when the patient was last monitored. Our approach is trained on batches comprising random image pairs I j , I k (Figure 1 ) of the same eye captured at different time points t j , t k ∈ R ≥0 from two visits, such that I j precedes I k with t j < t k . 
method: 
 AMD Stage Classifier: Both I j , I k are input to the same ConvNeXt-Tiny Encoder to obtain the features f j and f k respectively. ConvNext-Tiny with 29M parameters and 4.5GFLOPs is comparable to ResNet-50 and outperforms similar-sized Vision Transformer (ViT) architectures (Liu et al., 2022) , making it a suitable choice for our task. The stage classifier's GT y cls j = 1 if both T j ≤ 0 and E j = 1, otherwise y cls j = 0 for the scan I j . The stage classifier predicts the logit r j . The probability for the current AMD stage for I j is p j = σ (r j ), where σ(.) denotes the sigmoid activation. Notably, r j is proportional to the distance of f j from the decision hyperplane of the AMD stage classifier and would be used below to define a risk score for future conversion. The stage classifier treats each scan independently without considering any correlations between two scans of the same eye from different time-points. While it enables the learned feature to capture pathologies to distinguish dAMD from iAMD, it may fail to capture more subtle retinal changes indicative of how AMD will progress in the future (Appendix Figure 3(a) ).Time-Series Prediction: To address these issues, we incorporate a N-ODE based continuous time-series predictor called Stacked Multihead GRU-ODE (SMGRU-ODE) to model the future trajectory of AMD progression in the feature embedding using the current scan. SMGRU-ODE evolves the current feature f j over a (t k -t j ) time-interval to independently predict the future feature f k for time t k directly from the prior visit I j , while the actual feature f k is also obtained from I k . The encoder, SMGRU-ODE and stage classifier can now be jointly trained with the AMD stage classification task:EQUATIONwhere L bce (y, p) is the binary cross-entropy loss. p j , p k and pk are predictions from the stage classifier for the features f j , f k and f k respectively. The SMGRU-ODE architecture is detailed in Section 2.1 Intra-eye Consistency: For the disease progression trajectory predicted by the SMGRU-ODE to be consistent, the features f k and its stage prediction pk should match the corresponding f k and p k , obtained directly from I k . This consistency loss between the features (L cns-f tr ) and the stage predictions (L cns-rsk ) are defined as:EQUATIONThese losses combined with L cls ensure that the learned feature embedding (Appendix Figure 3 (a) vs (b)) not only characterizes the current disease stage but is also sensitive to the subtle retinal changes that capture the trajectory of the disease progression in the future. Since the future visit scans are unavailable at test time, only the top branch in Figure 1 highlighted in green is employed to obtain future predictions with the N-ODE. Risk Score Ranking: L cls ensures that the iAMD and nAMD samples lie on opposite sides of its decision hyperplane without imposing any ordering among iAMD cases. We envision a regularized feature manifold (Appendix Figure 3 (b) vs (c)) which correlates the risk of disease progression of a feature point to be inversely related to its distance from decision hyperplane, i.e., the closer an iAMD sample is to the decision hyperplane, the smaller its time to conversion, until it crosses over the hyperplane to the dAMD class. In this case, the logits r from the stage classifier acts as a risk score for AMD progression as it is proportional to the distance of the sample from the decision hyperplane. While predicting the probability of conversion within specified time-points (CDF) requires the N-ODE during inference, a scalar risk score is directly obtained from the logits of the stage classifier from the current scan. We consider a loss L rank which is defined using pairs of samples to encourage such ordering. Given a training batch comprising B pairs of images (i.e., a total of 2×B scans), we form all possible pairs (I m , I n ) which may or may not come from the same eye. L rank defines an auxiliary classification task where the difference of their scalar logits from the AMD stage classifier is fed through a neuron (with a single input and output) to obtain the probability of ranking r m > r n as P m>n = σ (w • (r m -r n ) + b), where w and b are scalar parameters of the neuron and the loss is defined as L rank = L bce y m>n rank , P m>n . The GT y m>n rank = 1, if T m < T n and E m = 1 (indicating the I m converts before I n ) or when T m < T n and both I m , I n are scans of the same eye (the risk increases in the future visits as damage to the retinal tissue is irreversible). Similarly, y m>n rank = 0 if T m > T n and, either E n = 1 or I m , I n come from the same eye, which signify cases where I n converts before I m . The image pairs that do not fall into either one of these two categories cannot be ranked due to censoring and are considered to have missing labels that are masked out during the loss computation. Finally, the total loss to train the proposed model is:EQUATIONwhere the loss weights λ 1 , λ 2 , λ 3 and λ 4 are not handcrafted but dynamically adapted during training using MTAdam (Malkiel and Wolf, 2021) (see Appendix E for more details). Upon projecting the initial scan I j to the feature f j , the N-ODE predicts its future trajectory in the feature embedding to model disease progression. Let f (t) represent the feature after a time t has elapsed since I j was imaged. As time progresses from t to t + dt by an infinitesimal amount, f (t) is displaced by v D • dt where v D denotes the instantaneous velocity vector. This can be modeled in continuous time using the N-ODE df (t) dt = v D (f (t)) with the initial value f (0) = f j , where v D is modeled with a DL network. We assume a time-invariant system, i.e., v D is solely dependent on the current feature f (t) and not on the time t elapsed so far.EQUATIONWe propose the SMGRU-ODE network to model v D (f ) which extends GRU-ODE (De Brouwer et al., 2019) (Dosovitskiy et al., 2020) and recent CNN architectures (Liu et al., 2022) , (Xie et al., 2017) . In Figure 2 (a), each layer employs f (t) as the hidden state and except for the first layer, also accepts an external input v d-1 (t) from its previous layer. Additive skip residual connections are applied between the inputs and outputs of each layer (v d (t) = vd (t) + v d-1 (t)). As depicted in Figure 2 (b), each head 1 ≤ h ≤ H independently projects the two, d-dimensional (768 for ConvNeXT-Tiny (Liu et al., 2022) ) inputs to a d/H dimensional sub-space using the fully connected (FC) layers, φ h,d (f (t)) and ϕ h,d (v d-1 (t)) which project f (t) and v d-1 (t) to f EQUATIONEQUATIONEQUATIONEQUATIONd (t) = 1 -u (h) d (t) ⊙ g (h) d (t) -f (h) d (t) (4b)whereΨ (rst)h,d (.) and Ψ (updt) h,d(.) comprise a FC layer followed by Layer Normalization (LN)and sigmoid activation to compute the update and reset gates, u The LN and activations are instead applied at the beginning of the next layer in φ h,d+1 (.) and ϕ h,d+1 (.). This is to ensure that (i) the output of the final layer v D can take arbitrary (including negative) values and (ii) the normalization and activation is applied after the additive residual connections in each layer (Figure 2(a) ) so that backpropagation gradients can be improved through pre-activation (He et al., 2016) . During forward pass of the N-ODE, the feature for a future time-point k is given by f k = f j + k-j 0 v D (f t+j )dt which can be numerically estimated using any black-box ODE-solver. During the backward pass, the computational graph related to each iteration of the ODE-solver is not saved but can be estimated on the run by solving another augmented ODE introduced by the adjoint sensitivity analysis in (Chen et al., 2018) . As a result, the training requires a constant amount of memory independent of the solver's step size and integration time allowing us to evolve the trajectory over long time-intervals, even with limited GPU memory. 
experiments and results: 
 Dataset: It consists of 3,534 OCT scans from 235 eyes (40 converters and 195 censored) from 123 patients, collected at the Department of Ophthalmology, Medical University of Vienna (Schlanitz et al., 2017) and acquired using a Spectralis scanner at a resolution of 49 B-scans (slices), each with a 512 -1024 × 496 pixels. Each eye was imaged every 3-6 months, with total follow-up periods spanning 2-7 years. For converter eyes, labels for each scan were computed by measuring the time interval between its acquisition and the first conversion visit. Our Pytorch code is available at https://github.com/arunava555/ Multihead_GRU_ODE_based_Survival_Analysis.Experimental Setup: A stratified five-fold cross-validation was performed by randomly dividing the scans at an eye-level to reduce the bias of a specific train-test data split. Each fold had 47 eyes with 8 converters and the number of scans varied between 667-707 across the folds. The model was trained five times, treating each fold as the test set, while the remaining dataset was randomly divided into 80% for training and 20% for validation. While the already converted dAMD scans were used during training, they were removed from the test set during evaluation. The performance was evaluated for predicting the conversion to dAMD within 6, 12, 18 and 24 months using the Area under the receiver operating characteristic curve (AUROC). Balanced Accuracy was used to assess binary predictions obtained by setting a threshold on the conversion probabilities at an optimal operating point determined from the validation set in each fold. Additionally, the Concordance Index (C-index) was used to evaluate the proposed risk score. It quantifies a model's ability to provide a reliable (inverse) ranking of the conversion time, taking censoring into account.Ablation Results: In Table 1 , we analyzed the effect of the depth and the number of heads in SMGRU-ODE. Either reducing the depth D from 3 to 1 (in row 1) while keeping H=12, or reducing H from 12 to 1, while keeping D=3 (in row 2) had an adverse impact on the performance at all time-points, both in terms of AUROC and Balanced Accuracy. The C-index also reduced from 0.777 to 0.744 in both cases. This justifies our incorporation of multiple layers and heads in GRU-ODE. From rows 3-6, we perform ablation on the loss terms. In row 3, we train the model with L cls loss (see eq. ( 1) ), which is the minimal loss required to predict future conversion without incorporating any other losses to regularize the feature embedding. Introducing L rank loss to it leads to a significant improvement in C-index (from 0.715 to 0.769) which is expected as L rank is geared towards improving the rank ordering. Moreover, it improves the conversion prediction performance for all timepoints both in terms of AUROC and Balanced Accuracy (except for AUROC -6 month). Next, introducing the L cns-rsk loss (in row 5) leads to further improvement in C-index, AUROC also improves for all except the 24-month time-point. However, the Balanced Accuracy shows mixed results with minor improvements for predicting conversion within 12 and 18 months but a slight drop in performance for the 6 and 24-month time-points. Finally, introducing the L cns-f tr loss leads to our proposed method in row 6. It consistently improves the AUROC, Balanced Accuracy and C-index metrics with the exception of the 6 month time-point. Overall, the results demonstrate the value of using all loss terms.Comparison with the State of the Art: In Table 2 , we compare our method against common survival analysis methods. 6-month time windows are considered for the discretetime survival models based on the censored cross-entropy loss (Wulczyn et al., 2020) and the logistic hazard model (Rivail et al., 2023 ). DeepSurv (Katzman et al., 2018) extends CoxPH with DL, while SODEN (Tang et al., 2022 ) is a N-ODE based method, previously used on tabular data. These methods were also trained with ConvNeXt-Tiny encoder but with modified classification layers and losses. Notably, all of these methods do not employ intra-  
conclusion: 
 A wide variability in progression speed and the lack of well-established biomarkers make predicting the progression of AMD challenging. We proposed a novel framework that combines an AMD stage classifier with a N-ODE to forecast dAMD onset at continuous future times. To learn meaningful features from scarce labels, we enforce (i) intra-subject consistency to ensure that the feature embedding is sensitive to temporal changes in the retina to predict the future; (ii) temporal ordering, where a scan's proximity to the AMD classifier's decision hyperplane is inversely related to its time-to-conversion. These constraints enabled our model to outperform several existing deep survival analysis methods. Additionally, temporal ranking allowed us to derive a scalar risk score to stratify eyes into low and high risk groups. While training uses longitudinal OCT scans, only a single scan at test time is needed for future conversion prediction. Our method for predicting dAMD onset can facilitate patient-specific disease management and enrich clinical trial populations with high-risk patients. Currently, the proposed method has been evaluated on a single-center dataset. Further evaluation of our method on multi-center data and adaptation to other survival analysis tasks in the medical domain, such as progression-free survival in cancer patients, are potential directions for future work. Use of a ViT based encoder and incorporating segmentations of relevant retinal layers and lesions as additional inputs may also be considered in the future to further improve performance.of 10 -6 and 10 -4 respectively. The performance on the validation set was monitored at the end of each epoch for early stopping with a patience of 50 epochs. The N-ODE was implemented with the torchdiffeq library (Chen, 2018). The Euler method was used as the ODE-solver due to its computational efficiency with a step size of 0.06, where the time between 0-3 years was mapped to [0,1]. During training, each batch was constructed with random training image-pairs (I j , I k ) with a time-interval of 0-3 years between them. The training batches were constructed to ensure that all I j were in the iAMD stage while half of the I k in each batch were in the dAMD stage (through oversampling) to enable the training of the AMD stage classifier. The proposed method required around 6 GB of GPU memory to train using a training batch size of 16 image pairs. The ConvNeXt-Tiny (Liu et al., 2022) encoder was initialized with the standard Image-Net pre-trained weights for end-to-end fine-tuning. The proposed SMGRU-ODE model with D=3, H=12 has 4,798,848 learnable network parameters.for each population group, computed across 1000 re-samplings using bootstrapping. The survival curves for the three risk groups show a clear separation, thereby demonstrating the effectiveness of the proposed risk score. The learned feature embedding (Fig. 4 (b) ) exhibit a smooth transition from fast(red) to slow converters(blue) along the feature manifold where the gray dots represent the censored scans.The Saliency maps obtained for the risk scores in Fig. 5 show the network to be sensitive to the structural changes around the RPE (e.g.  
 
 Appendix:  


Decision: This paper proposes a framework to predict AMD progression using neural ODEs. The paper was very well received and was recommended for oral presentation by all reviewers. Regrettably, the authors did not provide a rebuttal. Nevertheless, the reviewers mostly kept their original score and I continue to accept the paper for presentation at MIDL.

I strongly urge the authors to incorporate the reviewer comments in the camera-ready version and respond to reviewer comments in the future. I thank the reviewers for their valuable help.

Reviews: review: paper_summary: paper_summary: The paper explores disease progression using neural ODEs and makes a number of technical contribution paired with good experimental validation to predict iAMD to dAMD progression.
main_review: main_review: Since, they were a number of small concerns that I raised during the review, I am not very happy about not receiving a response. Nevertheless their seems to be an agreement that the paper should be accepted. So I hope the authors will incorporate some improvements in a final version
strength_weakness: strength_weakness: The relevance of censoring is mentioned in the beginning but not clearly resolved during the method and/or experimental section.
The result presentation is slightly underwhelming, the numerical tables have poor readability (very small font) and distribution plots or similar are missing (apart from Appendix Fig. 4). Some more qualitative results could be added there as well.
questions: questions: Some sentences/statements remain unclear: e.g. parametric models like CoxPH face limitations in hazard function distribution modelling flexibility. This is really hard to grasp for a non-expert Isn't the combination of parametric models with deep networks more capable of flexibel modelling? Furthermore, the concept of hybrid models (ones that extract segmentations from images) seems very useful but is not explored: could the authors envision further gains by incorporating this additional expert supervision?
paper_summary: paper_summary: The goal of this work is to develop an approach to predict the risk of age-related macular degeneration (AMD) from optical coherence tomography (OCT). They propose to combine an AMD stage classifier and a novel architecture using stacked multi-head GRU and a neural ODE to exploit intra-patient longitudinal information. OCT images. A siamese architecture trained on random image pairs from the same patient learn embeddings, a disease classifier and a prediction of the trajectory thanks to the N-ODE. The database consists in more than 3500 scans from one center. The validation of the approach through the ablation study and results (outperforming state of the art) are really convincing.
main_review: main_review: I would have appreciated a short revision to address my concerns on the redaction at least. Yet I stand by my initial rating on this work that deserves to be considered for acceptance : strong accept.
strength_weakness: strength_weakness: The work does not have any significant weaknesses:
- redaction could be improved (see detailed comments) to correct some small misprints and make it easier to read. 
- this is a single center study that would have to be confirmed in a multi-center study.
questions: questions: This would not change my mind, but some perspectives could be given in the conclusion.
paper_summary: paper_summary: The authors train a classifier to discriminate between intermediate and dry stages of Age-related macular degeneration (AMD), modeling the future trajectory of the disease progression in the learned embedding space. 
Interesting and well motivated constraints are imposed, including a temporal ordering so that the distance of a scan from the decision hyperplane of the AMD stage classifier is inversely related to its time-to-conversion. An intra-subject temporal consistency is also included. 
The proposed method is evaluated on a middle-size dataset (235 eyes) with AUC of predicting conversion within the next 2 years, and C-index. Extensive results, comparisons and analyses are provided.
main_review: main_review: I still think that minor improvements can be made, and that a response to all reviewers is important. However, my original recommendation is unchanged. I hope the authors will update the manuscript. Reviewing requires time (here 3 reviewers + AC), and it is important to incorporate this feedback even if the paper is likely accepted.
strength_weakness: strength_weakness: Nothing major. Figure 1 could be more self-explanatory, mentioning the main blocks and losses for a quick overview of the method (without making it too complicating).
No statistical tests are performed.
questions: questions: Nothing major, fix the minor comments and make the figure more self-explanatory.

Comments: title: comments to reviewers
Gentle reminder of rebuttal deadline
comment: Dear authors,

The PCs had access to all your answers, and they considered this for the final decision of the paper. Now, your answers are public to everyone, and no additional steps are required.

Thank you for your submission,
Maria, on behalf of the organization team
Dear Program Chairs,

During the rebuttal phase, I had selected the wrong settings by mistake and my replies were not visible to the reviewers. I realized it now after seeing the final response from the reviewers. I apologize for this mistake and I have now edited the settings so that everyone could see the replies. Please suggest  what additional steps could be taken now to address this mistake. 

Thanks and Best Regards,
Arunava Chakravarty
We thank the reviewer for recognizing the value of our contribution.  In the revised manuscript, we have made the following changes to address the reviewer’s comments.

1. We have now modified Fig. 1 as suggested by the reviewer to include the losses and extended the caption to make it more self-explanatory.


2. We agree with the reviewer that the proposed method should be evaluated on other survival tasks in the medical domain in the future.  We have now mentioned it at the end of the Conclusion section.


3. We have ensured that the abbreviations are defined at the place of their first use such as Age-related Macular Degeneration (AMD) in the Abstract  and Cumulative Distribution Function (CDF) in the introduction section. 


4. We have now also fixed the minor typos pointed out by the reviewer.

5. We have modified the abbreviation used for Concordance Index from CCI in the earlier version of the manuscript to C-index.
We are grateful to the reviewer for recognizing the merit of our paper. In the revised manuscript, we have made the following changes based on the reviewer's valuable suggestions.


1. We have expanded our conclusion with a discussion of some limitations of the current work which provide us with possible directions for extending the current work. This includes adapting our method on other survival analysis tasks in the medical domain and evaluating our current method on datasets acquired across multiple sites as pointed out by the reviewer.


2. The Balanced Accuracy was computed by selecting an optimal threshold on the validation set in each fold.  This has been clarified in the revised manuscript.


3. We have defined the full-form of the acronyms AMD and OCT both in the paper’s title  and Abstract as suggested by the reviewer to improve clarity. 


4. We have made minor modifications in the notations used in Fig 2 as suggested by the reviewer to improve readability.


5. We are grateful to the reviewer for pointing out the mistakes in the mathematical notations on pg. 5 and pg. 6 which we have now corrected in the revised manuscript.
We thank the reviewer for helpful suggestions and constructive feedback. Detailed discussion follows.

1. Vision Transformers (ViT) are known to require large datasets and suitable pre-training to avoid over-fitting. Moreover many ViT architectures are resource intensive with their computational complexity and memory requirements increasing quadratically with the number of input tokens making them unsuitable for our preliminary experiments with Neural-ODEs. However, we agree with the reviewer that Ablation experiments with different encoder architectures including ViTs need to be carried out as a part of future work which we have now mentioned at the end of the Conclusion Section. 


2. Reasoning behind using the ConvNeXt-Tiny architecture as backbone: ConvNext-Tiny with 29M network parameters and 4.5GFLOPs is comparable to ResNet-50 (25.6M parameters and 4.1GFLOPs). Hence, it is a reasonable choice for our medium sized dataset. Moreover, ConvNext-Tiny has been shown to outperform some of the ViT models of similar size on the Image-Net classification task such as Swin-Tiny, DeiT-Small and DeiT-Big. Compared to the ConvNext-Tiny Encoder, our proposed Multi-head GRU-ODE architecture is relatively small with 4.79M parameters. (in contrast to 29M for ConvNeXt-Tiny). We have now mentioned the motivation to use ConvNeXt-Tiny at the beginning of page 4 under the “AMD stage Classifier” heading. 


3. We thank the reviewer for pointing out a relevant work on survival analysis based on the Transformer architecture by Hu et. al. which we have now added to our related work. 


4.  “parametric models like CoxPH face limitations in hazard function distribution modelling flexibility” In DeepSurv and similar deep learning extensions of the CoxPH model, a single baseline cumulative hazard function $H_0(t)$ is computed for the entire training dataset using the Breslow Estimator. The cumulative hazard function for a specific eye with scan I, denoted as $H(t|I)$, is then determined as $H(t|I) = H_0(t) * \exp(r)$, where $r$ is a constant scaling factor across time. This means that $H_{0}(t)$ serves as a "prior" cumulative hazard across the population, scaled by $\exp(r)$ to generate individual cumulative hazard functions for each eye. Despite predicting $r$ with Deep Learning, the underlying assumption remains that hazard functions are scaled versions of $H_0(t)$. This constraint limits the flexibility of deep CoxPH models, as real-world hazard functions for different individuals may follow various monotonically increasing shapes, not necessarily scaled versions of a common $H_0(t)$.  In the revised manuscript we have now modified the line “parametric models like CoxPH face limitations in hazard function distribution modelling flexibility” to “Although CoxPH has been extended with DL using images \cite{katzman2018deepsurv}, they have not yet been explored to model AMD progression so far. Moreover, these models are inflexible as each patient's hazard function is constrained to be a scaled version of the same baseline hazard across the entire population.” 


5. “Furthermore, the concept of hybrid models (ones that extract segmentations from images) seems very useful but is not explored: could the authors envision further gains by incorporating this additional expert supervision?”     We agree with the reviewer's suggestion regarding potential performance improvement through the inclusion of additional inputs such as retinal layer segmentations and lesions, which we have now mentioned as part of future work in the conclusion section. However, this extension poses challenges, as it depends on the development of accurate segmentation networks for retinal layers and lesions in AMD patients. These networks may encounter segmentation errors and require collecting pixel-level segmentation Ground Truth for training. Additionally, the specific retinal layers and lesions that indicate future conversion risk remain clinically unclear, raising questions about which biomarkers should be prioritized. Furthermore, meaningful architectural adjustments are necessary in the encoder backbone to effectively process multiple inputs, including raw OCT scans and segmentation maps.


6. Additional qualitative results in the form of Gradcam Saliency Maps have been appended to the Appendix , Fig. 5.
Dear authors,

This is just a gentle reminder that the rebuttal deadline is coming up. I encourage you to take this opportunity so you can engage with the reviewers in the upcoming discussion phase.

Best wishes,
AC