Predicting Atrial Fibrillation Treatment Outcome with Siamese Multi-modal Fusion and Cardiac Digital Twins

Authors: alexander m zolotarev
abbas  khan
rayabat  khan
gregory  slabaugh
caroline h roney
Venue: MIDL_io_2024_Conference
Published: 2024-01-31 15:44:54

Abstract: Atrial fibrillation, the most common heart rhythm disorder, presents challenges in treatment due to difficulty pinpointing the patient-specific regions of abnormal electrical activity. 
While biophysical simulations of cardiac electrophysiology create a digital twin of atrial electrical activity based on CT or MRI scans, testing various treatment strategies on them is time-consuming and impractical on clinical timescales.
Our proposed pipeline, incorporating Siamese architecture, fuses latent representations of multi-modal features extracted from atrial digital twin before any therapy and predicts the outcomes of several treatment strategies.
A large in-silico dataset of 1000 virtual patients, generated from clinical data, was utilized to provide the biophysical simulations before (used for feature extraction) and after (used for calculating ground truth labels depending on whether atrial fibrillation terminates or not) various treatment strategies.  
By accurately predicting freedom from atrial fibrillation, our pipeline paves the way for personalized atrial fibrillation therapy with a fast and precise selection of optimal treatments.

Main paper: 
 introduction: 
 Atrial fibrillation (AF) remains the most common heart rhythm disease and is highly associated with an increased risk of stroke (Joglar et al., 2024) . Multiple clinical and experimental studies showed that AF may be caused and maintained by specific localized sources of pathological electrical activity (Hansen et al., 2015) . The locations of AF sources and their conduction pathways are highly dependent on the arrhythmogenic substrate, including fibrotic tissue distribution (Hansen et al., 2017) . Fibrosis is remodelled heart tissue with specific conduction parameters. It plays a crucial role in AF maintenance and can be used as a feature map to predict AF treatment outcome (Althoff et al., 2022) .The treatment of AF may be surgical and consists of ablation (physical destruction) of tissue areas with pathological sources of electrical activity. The baseline ablation strategy recommended by clinical guidelines is pulmonary vein isolation (PVI), which uses radiofrequency catheter ablation to electrically isolate the pulmonary veins from the rest of the heart, preventing triggers from the pulmonary veins from maintaining AF (Joglar et al., 2024) . However, PVI alone is usually not enough to successfully eliminate AF, and the long-term success rate for persistent AF ablation therapy is between 50 and 60%. There are multiple options for additional ablation lesions to stop the propagation of pathological electrical activity through the atria, including linear ablations, isolation of the left atrial appendage, or ablation of low-voltage areas and fibrosis (Rosa et al., 2021) .Human digital twins allow assessment of patient-specific features and simulation of physiological processes. They can provide personalised treatments incorporating information on both patient-specific anatomy and physiology. Specifically, the digital twin of the human heart can help to predict sudden cardiac death and to guide the optimal treatment for heart rhythm disorders (Trayanova et al., 2023 ) such as AF using biophysical simulations.However, AF biophysical models require a long simulation time, and predicting the longterm response to any therapy is challenging. Merging AF simulations with Deep Learning (DL) algorithms can overcome this limitation and allow us to predict the outcomes of various ablation strategies on clinical timescales. These predictions can base on the patient-specific characteristics of the electrical activity within the heart -by analyzing the atrial shape, its structure and parameters of electrical signals (anatomical and physiological features). As an example, fibrosis distributions are unique for each patient and affect the patterns of electrical activity. Late Gadolinium Enhancement (LGE) MRI is a non-invasive tool for estimating fibrosis quantitatively. Gadolinium accumulates in fibrotic tissue in higher amounts than in healthy tissue or the blood pool, and as a result, the intensity of LGE-MRI images from scar tissue is higher. Image intensity ratio (IIR, MRI intensity divided by the intensity of the blood pool) provides a method for quantifying atrial fibrosis distribution, standardized across different patients and MRI machines (Khurram et al., 2014) . Dominant Frequency (DF) and Phase Singularity (PS) maps have been proposed to be relevant attributes to characterize atrial electrical activity. Spectral analysis of atrial signals can locate the regions with the highest frequencies likely to maintain arrhythmia. In particular, DF maps visualize the frequency corresponding to the highest peak in the Fourier spectrum for each spatial location. Ablation of areas of high DF results in AF termination (Sanders et al., 2005) and the frequency features, including DF, may serve as successful indicators for localizing the area of AF by ML classifiers (Zolotarev et al., 2020) . As a further metric, the phase of transmembrane voltages can detect regions with rotational patterns which may sustain fibrillation. Specifically, phase singularities (which indicate the centres of rotational activity and areas of wavefront breakup) were shown to serve as a marker for the drivers which maintain AF pattern (Umapathy et al., 2010) .In this study, we propose a novel deep learning approach to predict the outcomes (the probability of atrial fibrillation termination) of different treatment strategies based on a large in-silico study of 1000 virtual patients. Our paper makes the following contributions:• We introduce a novel deep learning pipeline based on a Siamese architecture that fuses multi-modal feature maps derived from both atria.• We produce the methodology of generating multi-modal feature maps using an atrial digital twin before any treatment.• We show the effectiveness of the proposed pipeline and demonstrate that the full combination of multi-modal feature maps, Siamese architecture and fusion strategy is important to achieve the predictive capability. 
methods: 
 Dataset We used 1000 bi-atrial meshes from a statistical shape model derived from Cardiac CT scans of 19 healthy patients (Rodero et al., 2021) . In brief, right (RA) and left (LA) atrium were segmented from CT scans, and the surface meshes were utilized to form the bi-layer model of the heart. To create a more realistic representation of cardiac tissue, different anatomical structures were added to the mesh, such as inter-atrial connections to imitate the signal propagation between the atria, the crista terminalis and pectinate muscles, and the sinoatrial node, which is the heart's primary pacemaker (Roney et al., 2023) . We also mapped the fibre field to the atrial surfaces from the same bi-layer atlas to include heterogeneity within the tissue by highlighting the primary direction of signal propagation with the highest velocity and the perpendicular direction with a slower speed.To assign these structures, we utilised the universal atrial coordinates (UAC) (Roney et al., 2019) . UAC were calculated for each mesh by solving a Laplace equation with boundary conditions. This allows the comparison of spatial data across patients and the construction of a 2D representation of any 3D feature map. Next, each bi-atrial mesh was assigned two (for LA and RA) randomly selected fibrosis distributions from 100 clinical LGE-MRI scans of 100 AF (43 paroxysmal and 57 persistent) patients undergoing first ablation, as described in (Roney et al., 2022) . Based on the fibrosis distributions, the ionic and conductivity properties of fibrotic regions were changed to represent atrial remodelling. All of the aforementioned resulted in creating the in-silico dataset of 1000 virtual patients, which can serve as an initial foundation for biophysical simulations (Figure 1A ).The biophysical simulations were performed by solving the monodomain equation (Eq. 1, (Clayton et al., 2011) ) using openCARP software with the Courtemanche ionic model for the cellular action potential (Courtemanche et al., 1998) :EQUATION)where V is transmembrane potential, I ion -transmembrane ionic current per unit area, C m -the membrane capacitance per unit area, and D -diffusion tensor (m 2 s -1 ). AF was initiated by adding four Archimedean spiral waves on the atrial surface, with more details at (Roney et al., 2020) . AF episodes were simulated for 15 seconds to have enough information for phase and frequency analysis and feature extraction. Four types of ablation strategies were applied after 5 seconds of AF simulation across the virtual cohort, including PVI and PVI with LA, RA and bi-atrial fibrosis ablation. We applied ablation lesions by setting the tissue conductivity to 0.001 S/m for the corresponding elements of the atrial mesh (Bayer et al., 2016) . PVI ablations were simulated by adding two non-conducting rings around the left and right pulmonary vein antra.An automatic methodology was utilised to ablate the regions with higher fibrosis intensity. Firstly, IIR distributions for LGE maps were binarised with the threshold of 1.22 (Roney et al., 2022) . Regions of positive labels were connected to create several clusters of high fibrosis. Finally, each cluster was connected with non-conductive lesions to either PVI lesions or to heart valve borders (which are also non-conductive), depending on which was closer. The resulting lesion maze prevented the pathological electrical signals from freely propagating through the atrial surface and may lead to AF termination. Feature maps The DL pipeline was trained based on different input feature maps, including an ablation mask for each type of ablation and additional anatomical (fibrosis distribution), as well as physiological (Phase Singularity Density, PSD, and Dominant Frequency, DF) features calculated from the AF simulations before any ablation. AF simulations yield transmembrane voltage recordings with a temporal resolution of 50 ms for all nodes on the atrial mesh. The DF values for each mesh element were calculated as the frequency of the highest spectral peak of the signal, excluding peaks of the spectrum above 20 Hz to remove physiologically uninformative frequencies from the analysis. The phase of atrial signals was calculated by applying the Hilbert transform after subtracting the signal mean. We constructed PSD maps by counting the number of PS occurrences (by calculating the topological charge (Rogers, 2004 )) across the whole duration of the AF simulation at each mesh element and smoothing the values through the mesh.For both PSD and DF values on a 3D atrial mesh, we used the UAC to convert them from 3D to 2D feature maps. Initially, we set up the shape of 2D maps, which was selected as 128 by 128 pixels, with a normalized coordinate system (0-1) on both axes. We calculated the 2D UAC coordinates for each mesh element on LA and RA surfaces, where both the x and y axes also span between 0 and 1. Each pixel of the 2D map corresponds with a mesh element with a minimal distance between the pixel centre and the 2D location of this mesh element. By assigning DF and PSD values from selected mesh elements to corresponding 2D pixels for both LA and RA, we generated four 2D feature maps: DF LA, DF RA, PSD LA, PSD RA. Finally, we created 2D ablation masks and fibrosis maps based on the same protocol as for DF and PSD maps (Figure 1A ). Labels of ablation lesion elements from the 3D atrial mesh served as the basis for ablation masks (mask LA, mask RA), while fibrosis distributions from LGE-MRI scans were used to create fibrosis feature maps (F LA, F RA).The binary outcomes of simulations after different ablation strategies were used for training the DL pipeline and were assigned a ground truth (GT) label, whether AF terminates or not. We counted the AF episode as terminated if there is no electrical activity at the end of the recording (specifically, the last peak of the action potential is within the range 0 -60% of the whole duration of the recording). To define the final GT, we calculated the last peak's values of the AF signal for each node of both LA and RA and averaged the value across all nodes (Figure 1B ).To address the problem of fast and correct prediction of AF ablation, we propose a DL pipeline based on Siamese architecture (Figure 2 ). Each head consists of DenseNet121 (Huang et al., 2016) network and utilises 4 channel-wisely concatenated feature maps (ablation mask, PSD, DF and fibrosis maps, all of the size 128 by 128 pixels) from LA and RA separately. The outputs of both heads (n = 32 each) were fed into a Multi-modal Outer Arithmetic Block (MOAB) (Alwazzan et al., 2023) for the fusion of latent representations of features. It allows us to capture and combine the relevant features from different anatomical structures by applying four arithmetic operations. The resulting arrays were concatenated into multi-modal tensor which passed through the 2D convolution layer. The final prediction of AF ablation outcome as the probability of AF termination was achieved through two successive fully-connected layers (n = 1089 and 512 respectively) and one dropout layer.The performance was evaluated by accuracy, ROC-AUC and F1-score on a testing set of 250 virtual patients. To avoid data leakage, we assigned the virtual patients with the same fibrosis distribution randomly selected from 100 LGE-MRI maps to either the training or testing set. The training was conducted on an NVIDIA GeForce RTX 3080 video card using Binary Cross Entropy loss and Adam optimizer with a learning rate of 4e-7.  
results: 
 We show the effectiveness of the proposed model on the dataset of 1000 virtual patients described in Section 2. For each virtual patient, we ran the biophysical simulations with AF initialization before and after 4 types of ablation strategies. All 1000 AF simulations before ablation were sustained (AF was stable during the whole duration of recording). Overall, the percentage of successful ablations (which lead to AF termination) was 41%. Analysing this by ablation therapy approach, the effectiveness of PVI ablation was 25%, PVI with LA fibrosis ablation -36%, PVI with RA fibrosis ablation -39%, and PVI with bi-atrial fibrosis ablation -60%. An example of one virtual patient with only PVI with bi-atrial fibrosis ablation resulting in AF termination is shown in Figure 3 . We demonstrated that the proposed model achieved high performance in predicting the outcome of different ablation strategies, including a ROC-AUC of 0.92, an accuracy of 87%, and an F1-score of 0.88. In particular, the outcomes of PVI ablation alone are the easiest to predict, with ROC-AUC of 0.93 and 21 incorrectly predicted cases (out of 250). In contrast, 44 cases of PVI with RA fibrosis ablation were predicted incorrectly, which resulted in a ROC-AUC of 0.904 for this type of ablation. Specifically, the confusion matrices for predictions in general and for each ablation type are provided in Appendix A. The running time for one ablation strategy prediction was 41 seconds (34 -for the feature extraction and 7 -for the inference), which means that the DL pipeline can predict AF treatment outcome more than 3 times faster in comparison with conventional biophysical simulation (the fastest model run was at least 134 seconds on the UK national supercomputer ARCHER2).We conducted an extensive ablation study to test the individual effects of each component of the proposed DL pipeline (Table 1 ). Firstly, we checked whether it is necessary to use two heads by comparing it to the model's performance if we used a one-headed pipeline (all feature maps (n=8) were initially concatenated and fed into the DenseNet121 network to provide the final prediction after passing fully-connected layer and softmax activation). Then, we tested what would happen if we changed only one feature from the best implementation (two DenseNet121 heads with MOAB fusion of outputs). We utilised the different backbone architecture for two heads (EfficientNet) rather than DenseNet121, then the simple concatenation of outputs rather than MOAB fusion, and have shown that the original model has the best metrics. Finally, we showed that separating feature maps before training into LA and RA parts (which is physiologically relevant) performs better than separating them into 2 ablation masks vs all (fibrosis, PSD and DF maps) feature maps.We tried to answer the question of how different input feature maps affect the performance of the proposed model. For comparison, we selected the best-performing implementation consisting of two DenseNet121 heads with input features from the left and right atrium and outputs fused by MOAB. We found that the baseline model trained solely on the ablation mask acquired a ROC-AUC of 0.69, whereas incorporating fibrosis, DF or PSD elevated the metric to 0.86, 0.87, and 0.91, respectively. Furthermore, combining all features led to a ROC-AUC of 0.92 ( Discussion A recent multi-centre clinical study showed that ablating fibrosis in left atrium (LA) after PVI does not change the rate of AF recurrence in comparison with PVI alone (Marrouche et al., 2022) . However, this study investigated the effect of fibrosis ablation targeting the LA only, while it has been shown by an experimental study of adenosine-related AF (Li et al., 2016) that fibrotic tissue in right atrium (RA) can also serve as an anchor for AF pathways. Therefore, RA fibrosis areas could serve as an ablation target, and we have shown that ablating these areas improves the outcomes in a large insilico trial. We have previously demonstrated that machine learning classifiers can predict long-term AF recurrence rates for 100 AF patients with high accuracy. Also, we found that ROC-AUC values varied from 0.61 if relying only on patient history to 0.85 if adding imaging and simulated physiological inputs (Roney et al., 2022) . Confirming these findings we have currently shown that selecting input features is crucial and adding physiological characteristics such as PSD improves the results.Limitations Our study is not without limitations. Firstly, all findings are based on synthetic data and further validation on clinical recordings is needed. Secondly, all atrial meshes were mapped by the same fibre field and the effect of adding different fibre maps on the same atrial mesh can be our future task. The area near the sinoatrial node was not excluded from the RA ablation masks. 1000 meshes were generated from a limited number of atrial anatomies, and although it was shown that the anatomical volumes of this virtual cohort were within the range of values from UK Biobank cohort of over 5000 individuals (Petersen et al., 2017) , the generalizability issue should be tested more carefully on a bigger set of clinical meshes. The threshold for fibrosis distributions and the resolution of feature maps were set up to be 1.22 and 128 respectively, and changing the values of these hyperparameters may affect the results. The current simulation data used for training our pipeline assumes homogeneous electrophysiological properties in tissue conductivities used across healthy or fibrotic tissue. However, fibrotic remodelling is multifactorial, which may affect AF properties (Roney et al., 2016) . Lastly, fibrosis distributions were randomly assigned to each of the LA and RA separately meaning that virtual patients may have a LA fibrosis distribution from a paroxysmal patient and RA from a persistent one, or any such combination. Our future work will investigate the effects of AF classification on fibrosis distribution and simulated AF as well as the effects of type of remodelling on prediction. 
conclusion: 
 This study presents a novel approach for predicting the outcome of atrial fibrillation ablation by integrating an atrial digital twin with a deep learning pipeline. The proposed pipeline leverages a Siamese architecture with multi-modal fusion of latent representations, resulting in superior performance. Furthermore, by incorporating anatomical and physiological features from pre-ablation atrial fibrillation simulations into the model, we demonstrate significant improvement in predictive capability compared to solely using an ablation mask.  
 
 Appendix:  


Decision: This paper received the ratings with 1 borderline, 1 weak accept and 1 strong accept. The authors have addressed the major concerns during the rebuttal. Therefore a decision of accept is recommended. I strongly suggest the authors take the reviewers’ comments into the final revision.

Reviews: review: paper_summary: paper_summary: The authors present a method to assess the success of treatment ablation based on feature maps coming from imaging data and biophysical simulations using a digital twin created from CT or MRI scans. The long-term response of the treatment is predicted using a siamese neural network without the need of running long simulation times.
main_review: main_review: The authors have clarified all my doubts, and my concerns were properly addressed. Furthermore, they did a great effort updating the paper and providing further evidence as well as detailed comments. I think it is a good paper, and I am confident it will be of interest for the community.
strength_weakness: strength_weakness: - In this in-silico dataset the success rates of the tested ablation procedures are in the range [25%-60%]. Thus, the best procedure including LA and RA fibrosis maps achieves the success rate of PVI in clinical practice. Therefore, it seems that it could be easier to bias the network to predict mostly not successful ablation, specially in the PVI case for the in-silico data. Would it make sense to balance the training between successful and not successful cases? Would it be possible to include a confusion matrix for the predictions?
questions: questions: - The fiber directions are mapped from a reference atlas. Thus, I assume that most of the subjects present a similar distribution. What is the effect of fibers direction? , I also assume that they play an important role in the conduction dynamics. From Table 2, PSD seems to be the most relevant feature map for the prediction, while fibrosis not so much, might it be related to the fibers direction? i.e: maybe the network is learning which 'fibers' ablate to better stop the abnormal patterns.

- In Section 3.1 the authors state that they sample 2 fibrosis distributions, on for each atrium, do these distributions come from the same sampled subject?
paper_summary: paper_summary: This paper presents a Siamese network for fusing multiple modalities of data, including those output from digital simulations, to predict the outcome to 4 different A-fib treatment. The input to the model includes PSD, DF, and fibrosis maps along with the ablation mask representing ablation strategies. The output is a binary prediction of whether the ablation is successful. Experiments are conducted on virtual patients, showing the ability of the model to predict outcome in held-out patients.
main_review: main_review: The authors addressed some of my questions. Given the challenging nature of clinical validation for this type of works, and the extensive amount of works put forth in the current work, I'm happy to recommend the acceptance of this work. It could raise good discussion within the MIDL community regarding how digital twins and AI can benefit each other.
strength_weakness: strength_weakness: The lack of any indication of the clinical validity of the prediction is a main weakness of the present study. 

The sample size in terms of the unique patients contributing to the basis of the virtual patient cohort is somewhat limited — 19 healthy patients. The 1000 bi-atrial meshes are obtained by combining these meshes with 100 different fibrosis distributions. How would this limit the variability of the simulation and thus the generalizability of the subsequent outcome prediction model is not clear.
questions: questions: It’d be helpful if the authors can respond to the questions about the possible clinical validation no matter how limited, as well as how the split of virtual patients is determined based on the atrial mesh used; comments are generalizability due to the limited unique atrial geometry and fibrosis samples should also be added.
paper_summary: paper_summary: The authors present a novel framework to predict the outcome of ablation therapy to treat atrial fibrillation (AF). Specifically, a digital twin of the left and right atrium is incorporated to provide electrophysiological features and training data for the deep-learning-based classifier. For each the left and right atrium, the dominant frequency and phase singularity map are derived from AF simulations before ablation and combined with a fibrosis and an ablation map to form the features for a Siamese network build on DenseNet. The proposed pipeline was evaluated in two ablation studies by training it on 750 geometries and evaluating it on 250 unseen geometries. On this synthetic dataset, the proposed method demonstrates high accuracy in classifying whether the AF terminates under a specific ablation pattern.
main_review: main_review: The authors have clarified all questions and concerns in great detail (big thank you to the authors!). The changes made to the manuscript are appropriate and further boost the quality of the already good manuscript. In addition to general improvements in text, the proper discussion of the results as well as the additional insights through the confusion matrices greatly improve the paper. I am confident that this submission will be of great interest to the community.
strength_weakness: strength_weakness: - The proposed pipeline is only trained and evaluated on synthetic data 

- The training data seems to include only simulations with fixed electrophysiological parameters (properties of the ionic model and conductivity of fibrotic and healthy tissue), which may not reflect the inter-patient tissue variability 

- In both the method and results section, the paper misses several details to fully understand and assess the contribution (see detailed comments)

- The paper lacks a critical discussion of the results
questions: questions: In addition to the points raised under weaknesses and in the detailed comments, could the authors please comment on what role the resolution of the raw data, the mesh, and the feature maps might play onto the results?

Comments: title: Follow-up question
comment: Thanks for the clarification. 

Although more extensive and clinical validation are required in the future works, I appreciate the extensive amount of works involved in this proof of concept and acknowledge the difficulty in clinical validation. I'm happy to recommend the acceptance of this MIDL paper and look forward to seeing the future iterations of this interesting work.
Hi, if I figure out the question correctly,  
1000 synthetic geometries were split into 75% training and 25% testing sets, which means that 250 virtual patients in the testing set are based on 250 unique SSM geometries. All virtual patients are bi-atrial models based on 1000 SSM and 100 LGE distributions from LA. We did not use these real LA geometries anywhere, because we do not have RA geometries for them to construct bi-atrial models.
Thanks for the clarification. 

Outside the issue of clinical validation, can the authors clarify if the training/test split considered leaving out some LA geometry? It was not clear in the response above.
(#3.13) In addition to the points raised under weaknesses and in the detailed comments, could the authors please comment on what role the resolution of the raw data, the mesh, and the feature maps might play onto the results? 

We thank the reviewer for this valuable comment. The resolution significantly affects the results. Specifically, we did an additional ablation study to test how model performance depends on the resolution of the feature maps. We found that ROC-AUC scores were 0.908, 0.903, and 0.917 for resolutions of 64, 96, and 128 pixels. We did not test the effect of mesh resolution on the performance on 1000 SSM but another study from our group has shown that downsampling of the mesh used for analysis from the high-resolution simulation mesh can lead to a large change in phase singularities distributions (https://www.touchcardio.com/imaging/journal-articles/27-spatial-resolution-requirements-for-analysing-atrial-fibrillation-wavefront-patterns-insights-from-virtual-patient-models/). However, further investigations are needed to determine the optimal balance between data resolution, prediction metrics and run times. We have added this information to the Limitation section: “The threshold for fibrosis distributions and the resolution of feature maps were set up to be 1.22 and 128 respectively, and changing the values of these hyperparameters may affect the results.”).
(#3.9) Evaluation: Since the effectiveness of the ablation is varying significantly between PVI only and PVI + biatrial fibrosis ablation, it would be very valuable to understand to what degree this imbalance impacts the results, e.g., by showing confusion matrices for each ablation type. 

We thank the reviewer for this comment and want to highlight that we have not currently separated the training set into 4 datasets of different ablation strategies (all of them are unbalanced). The training set has a combination of all ablation strategies and has a good balance between positive and negative values (41:59). We found that confusion matrices for each ablation strategy will be helpful for the readers and have added them to the Appendix (as well as the general confusion matrix with respect to Q#1.1). 

(#3.10) One argument for the proposed pipeline is a significantly faster selection of the optimal ablation strategy compared to running several AF simulations with varying ablation protocols. Could the authors please comment on the run-time of the proposed method? How long does the inference of the DL network take in comparison to running the "baseline" AF simulation, i.e., before ablation, and extracting the electrophysiological features from the simulation?  

We thank the reviewer for pointing out this important information. To answer this question, we compare the running times of biophysical simulation and DL prediction for one hold-out case, which was PVI with bi-atrial fibrosis ablation for one virtual patient. The running time of AF ablation simulations was 351 seconds while running on the local computer with 64 core desktop CPUs and 134 seconds on the UK national supercomputer ARCHER2. In contrast, the running time for Deep Learning prediction was 41 seconds (34 seconds for extracting the feature maps and 7 seconds for the inference on NVIDIA GeForce RTX 3080 video card). Therefore, our proposed model can predict the ablation outcome from 3 to 8 times faster than conventional biophysical simulation, depending on the power of the CPU for simulations.  We have added this information to the results: “The running time for one ablation strategy prediction was 41 seconds (34 seconds for the feature extraction and 7 seconds for the inference), which means that the DL pipeline can predict AF treatment outcome more than 3 times faster in comparison with conventional biophysical simulations (the fastest model run was at least 134 seconds on the UK national supercomputer ARCHER2).”. 

(#3.11) Table 1: Could the authors please provide more information on the ablation experiments? Was special care taken so that each method has approximately the same number of parameters and an adequate latent dimension size before the fully-connected prediction head? On the experiment to separate the ablation masks vs all other masks: Was this network also trained as a Siamese network? If it was trained as a Siamese network, how was the input dimension of the 2 ablation masks adapted to match the 6-dimensional input of the other head? In addition, how can this be justified since I would assume that each head must learn significantly different features? 

Thank you for this interesting question. We have not currently tested any modifications of the network architecture related to ablation studies. During ablation masks vs all feature maps the network was trained in a Siamese (but not symmetric) mode, with 2 input channels for masks and 6 channels for feature maps (in contrast with 4 channels for the left and 4 channels for the right atrium as a default). The number of output features after passing the DenseNet121 remained the same as in the original MOAB paper (n=32) - for both fusion between 4 and 4 channels and the fusion between 2 and 6 channels. 

(#3.12) Table 2: The results of using PSD + mask are very close to using all features. Could the authors please comment on whether the improvement is statistically significant? 

To respond to this question, we compare the metrics on 4-fold cross-validation between the full set of input feature maps and PSD maps and ablation masks only. The resulting ROC-AUCs are shown in the Table below. The student t-test indicates that there is no statistically significant difference between these two sets of metrics (p-value = 0.08).  

full  vs PSD + mask:

1st fold - 0.917 vs 0.914
2nd fold - 0.904 vs 0.876
3rd fold - 0.933 vs 0.907
4th fold - 0.863 vs 0.856
(#3.1) The proposed pipeline is only trained and evaluated on synthetic data  

We thank the reviewer for this comment. By using synthetic data, we have more control over the impact of different feature maps including fibrosis distributions, we were able to test different ablation approaches, and it enabled us to choose and validate the optimal configuration of a Multi-modal Fusion approach. At a later stage we wish to validate our approach on clinical data.  

(#3.2) The training data seems to include only simulations with fixed electrophysiological parameters (properties of the ionic model and conductivity of fibrotic and healthy tissue), which may not reflect the inter-patient tissue variability  

We thank the reviewer for this question. In this study, we utilised different distributions of fibrosis to capture different distributions of conduction velocity and effective refractory period in our simulations. Our future work will investigate the effects of different types of fibrotic remodelling and restitution heterogeneity on our findings and whether the addition of these heterogeneities require additional input features. This will open new research directions such as what spatial resolution of restitution data is required for prediction. We have added the following to the limitations of our paper: “The current simulation data used for training our pipeline assumes homogeneous electrophysiological properties in tissue conductivities used across healthy or fibrotic tissue. However, fibrotic remodelling is multifactorial, which may affect AF properties (https://pubmed.ncbi.nlm.nih.gov/28011842/). Our future work will investigate the effects of AF classification on fibrosis distribution and simulated AF as well as the effects of type of remodelling on prediction.” 

(#3.3) In both the method and results section, the paper misses several details to fully understand and assess the contribution (see detailed comments) 

Thank you for noticing these omissions, we have replied to each of comments below. 

(#3.4) The paper lacks a critical discussion of the results 

We thank the reviewer for this comment and have added the separate discussion section to our manuscript. 

(#3.5) Introduction: While all statements in the introduction are correct, I believe that it would greatly help newcomers to the field to add few more citations that support the statements, e.g., please add to the statement "PVI is the baseline ablation strategy recommended by clinicial guidelines" a citation to one of these guidelines. 

We agree with the reviewer and add the reference to the most recent guideline, Joglar 2023, section 8.4.2. 

(#3.6) Equation 1: Shouldn't the diffusion term be divided by C_m? If I recall it correctly, the monodomain equation is also missing the variable to denote the membrane surface-to-volume ratio. Please double-check the equation. 

We used the form proposed in Clayton 2010 paper, equation 10 (https://doi.org/10.1016/j.pbiomolbio.2010.05.008). For more clarification, we have added this reference to the sentence with the equation and the units of diffusion tensor (m^(2) s^(-1)). 

(#3.7) Section 3.1, last paragraph: Could the authors please comment on whether such a lesion maze is clinically plausible? 

Yes, the ablation lesion set for either PVI or PVI plus any fibrosis ablation is clinically relevant. Specifically, we refer to the DECAAF II clinical trial, where the operator either encircled or covered with ablation lesions all fibrotic areas (10.1001/jama.2022.8831, section Fibrosis-Guided Ablation). However, we need to state that these ablation lesion sets cannot be directly applied to clinical patients. For example, we did not consider the ease of ablating different regions and we did not exclude lesions occurring in the vicinity of the sinoatrial node, the primary pacemaker of the heart, even though a clinician would be careful to not ablate these regions. We have added this information to the Limitations section: “The area near the sinoatrial node was not excluded from the RA ablation masks.”. 

(#3.8) Training / test data: Could the authors please confirm (or clarify) that the 1000 synthetic geometries were split into 75% training and 25% testing? Was an additional validation set used to fine-tune the hyperparameters? How were the shapes sampled from the shape model? Were any special considerations taken to sample the test shapes? 

Yes, we confirm that 1000 synthetic geometries were split into 75% training and 25% testing set, we didn’t fine-tune the hyperparameters on the additional validation set. The cases for the testing set were selected following the rule that a given fibrosis distribution can be within either the training or testing set.
(#2.1) Understandably clinical validation on this type of setup is difficult — nevertheless, on the actual patients that contribute atrial mesh or fibrosis map to the virtual cohort, there must be history of ablation treatment and outcome. Can this be used as retrospective validation of the proposed prediction network, at least as a proof of concept. 

The subjects that were used for constructing the anatomical statistical shape model do not have any history of AF treatment as they are healthy asymptomatic patients who went to the emergency with acute chest pain and had no cardiac conditions during follow-up. To incorporate fibrosis in the models, we used the LGE data of 100 patients undergoing ablation therapy for AF from a previous study (10.1161/CIRCEP.121.010253) where only left atrial simulations were performed. In our current study, we decided to use bi-atrial simulations to include the effects of the right atrium and right atrial fibrosis. We are currently working on constructing bi-atrial versions of these 100 cases for clinical validation. Please see more discussion about it in the response to Q#2.3 where we also point to new text in the manuscript. 

(#2.2) While the virtual patient is left out based on their fibrosis distribution, what about the atrial mesh? 

Please see our response to the previous question. Briefly, we had only left atrial (LA) segmentations and could perform AF biophysical simulations only on LA meshes; however, we are working on constructing bi-atrial versions of these 100 cases for clinical validation. 

(#2.3) The sample size in terms of the unique patients contributing to the basis of the virtual patient cohort is somewhat limited — 19 healthy patients. The 1000 bi-atrial meshes are obtained by combining these meshes with 100 different fibrosis distributions. How would this limit the variability of the simulation and thus the generalizability of the subsequent outcome prediction model is not clear. It’d be helpful if the authors can respond to the questions about the possible clinical validation no matter how limited, as well as how the split of virtual patients is determined based on the atrial mesh used; comments are generalizability due to the limited unique atrial geometry and fibrosis samples should also be added. 

The reviewer indicates a very important aspect of our work on how to translate these findings into clinical practice or at least show its clinical validity. Our ultimate goal is clinical validation of the proposed pipeline, and we present in this paper the first step of developing our approach on simulated data. We hope that the pipeline even without clinical validation can be interesting for data science and cardiac modelling communities and can serve as a proof-of-concept example of merging a novel deep learning architecture and a large in-silico trial of cardiac digital twins with different AF treatment strategies. By using synthetic data, we have more control over the impact of different feature maps including fibrosis distributions, we were able to test different ablation approaches, and it enabled us to choose and validate the optimal configuration of a Multi-modal Fusion approach. 

We are currently working on extensive clinical validation by adding a separate hold-out testing set of clinical patients as well as creating a more representative second version of the Statistical Shape Models based on additional CT/MRI scans for cases with AF. However, this process is time-consuming and requires significant development; as such, we unfortunately cannot perform new AF simulations for testing now. Once we have completed the bi-atrial segmentation of these cases with known ablation approaches and outcomes, we will test our deep learning pipeline predictions against the known outcome.  

Focusing on the generalizability issue, we want to emphasize that the authors of the original paper with 1000 Statistical Shape Models have shown that the anatomical volumes of the CT cohort were within the range of values from the UK Bio-bank in a population cohort of over 5000 individuals [https://doi.org/10.1186/s12968-017-0327-9 ]. Regarding the concerns of clinical validation and generalizability, we have added the following sentences to limitations section: “Firstly, all findings are based on synthetic data and further validation on clinical recordings is needed. 1000 meshes were generated from a limited number of atrial anatomies, and although it was shown that the anatomical volumes of this virtual cohort were within the range of values from UK Biobank in a population cohort of over 5000 individuals, the generalizability issue should be tested more carefully using a bigger set of clinical meshes.”.
(#1.5) The fiber directions are mapped from a reference atlas. Thus, I assume that most of the subjects present a similar distribution. What is the effect of fibers direction? I also assume that they play an important role in the conduction dynamics. From Table 2, PSD seems to be the most relevant feature map for the prediction, while fibrosis not so much, might it be related to the fibers direction? i.e: maybe the network is learning which 'fibers' ablate to better stop the abnormal patterns. 

We have mapped the same fibre field from a bi-atrial bilayer atlas for all virtual patients; therefore, with the current set of simulations, fibre orientation cannot be used as a valuable feature for learning. However, fibre orientation plays a role in signal propagation and the effect of adding different fibre maps on the same atrial meshes can be our future task. To further clarify this point, we changed the sentence in the Methods section (“We also mapped the fibre field to the atrial surfaces from the same bi-layer atlas to include heterogeneity within the tissue by highlighting the primary direction of signal propagation with the highest velocity and the perpendicular direction with a slower speed.”) and added the information about it to Limitations (“All atrial meshes included the same fibre field and the effect of adding different fibre maps will be investigated in future studies .”). 

(#1.6) In Section 3.1 the authors state that they sample 2 fibrosis distributions, on for each atrium, do these distributions come from the same sampled subject? 

These distributions were randomly assigned to atrial meshes and can be either from the same or different patients. We have stated this in the Methods section: “Next, each bi-atrial mesh was assigned two (for LA and RA) randomly selected fibrosis distributions from 100 clinical LGE-MRI scans”. We checked the virtual cohort and found that 1% of all patients have fibrosis distributions on LA and RA which came from the same patient.
R#1 

(#1.1) In this in-silico dataset the success rates of the tested ablation procedures are in the range [25%-60%]. Thus, the best procedure including LA and RA fibrosis maps achieves the success rate of PVI in clinical practice. Therefore, it seems that it could be easier to bias the network to predict mostly not successful ablation, specially in the PVI case for the in-silico data. Would it make sense to balance the training between successful and not successful cases? Would it be possible to include a confusion matrix for the predictions? 

We thank the reviewer for this valuable comment. The proportion of positive and negative values in the overall dataset is 41:59 (first paragraph of Results), which means that our dataset is well-balanced. We would like to indicate that we trained the model on all ablation options simultaneously and did not aim to separate datasets based on the ablation type. We also agree that adding the stratification into the dataloader for training is a good idea in the case of unbalanced datasets. This will become important if we look at the different ablation approaches separately (e.g. in the case of PVI ablation alone). We also decided to add to the Appendix a table with the confusion matrix for the prediction and link to this Table to the main manuscript. 

(#1.2) What’s the balance between paroxismal and persistent maps assigned to the in-silico data?  

The balance between fibrosis distributions from paroxysmal and persistent AF patients is 43 to 57 cases, as stated in section 3.1. However, we did not create the virtual patients based on only paroxysmal or persistent maps, these distributions were assigned to each of the LA and RA randomly to generate virtual patients with a range of fibrosis distributions. As such, the virtual patients may have a fibrosis distribution in the LA from a paroxysmal case, and in the RA from a persistent case, or indeed any combination. This choice means that the models capture different degrees and distributions of fibrosis observed across the spectrum of AF patients, but we cannot label virtual patients as paroxysmal or persistent. The effect of the AF classification (paroxysmal or persistent) on fibrosis distributions and resulting AF simulations is an interesting question and can be a task for investigation in future research. We have added the following to our limitations and future work: "Lastly, fibrosis distributions were randomly assigned to each of the LA and RA separately meaning that virtual patients may have a LA fibrosis distribution from a paroxysmal patient and RA from a persistent one, or any such combination. Our future work will investigate the effects of AF classification on fibrosis distribution and simulated AF.”  

(#1.3) Future steps or limitations of the current work in the conclusion would be appreciated. 

We appreciate this comment and have added the relevant information to the manuscript section of limitations. 

(#1.4) The authors used a fixed threshold for the fibrosis maps. Nevertheless, to my understanding, the fibrosis quantification from LGE is semi-quantitative, i.e: the values are scaled with respect to healthy tissue. Thus, subjects with high levels of fibrosis would not have a good 'healthy' reference and the amount of fibrotic tissue might not be estimated correctly. If that's the case, how could this threshold impact the results? 

We thank the reviewer for this question. Fibrosis intensities were normalized by the mean blood pool intensity to account for differences in image contrast. This scaling is more robust and less subjective than normalizing with respect to healthy tissue. However, the selection of a threshold value for the regions in which fibrosis was incorporated in the models has a significant impact on the results. To provide more clarity we have stated the information about normalization in the manuscript (“Image intensity ratio (IIR, MRI intensity divided by the intensity of the blood pool https://doi.org/10.1016/j.hrthm.2013.10.007) provides a method for quantifying atrial fibrosis distribution, standardized across different patients and MRI machines.”). We also added the sentence about the threshold to the Limitations section (“The threshold for fibrosis distributions and the resolution of feature maps were set up to be 1.22 and 128 respectively, and changing the values of these hyperparameters may affect the results.”).
We thank the editor and the reviewers for their comments and questions. We have provided comprehensive answers and have changed the manuscript to address their points; specifically, we removed the Abbreviations section and added new Discussion, and Limitations sections and Appendix. Regarding the strictly limited size of the manuscript, we have also slightly changed the text; however, all findings, tables and figures remain the same.