Personalized Evaluation of Atrial Complexity of Patients Undergoing Atrial Fibrillation Ablation: A Clinical Computational Study

Simple Summary Atrial fibrillation is a type of arrhythmia that occurs when the electrical activity of the heart in the atrium is not coordinated, and its consequences can be lethal. The driving source that initiates this chaotic activity can be located anywhere in the atrium, but most frequently appears in certain areas such as the pulmonary veins. In this study, we developed a new estimation method to evaluate possible source location and complexity of the arrhythmia using computer simulations. This method represents mathematical descriptions of natural processes that can be used to mimic a real scenario, including specific information such as the atrial anatomy. Here, we identified a specific biomarker the enabled obtaining a foci distribution map and found that elimination of pulmonary vein drivers was associated with a successful long-term ablation outcome. This study could, therefore, help to identify and characterize patients in order to better plan the ablation procedure. Abstract Current clinical guidelines establish Pulmonary Vein (PV) isolation as the indicated treatment for Atrial Fibrillation (AF). However, AF can also be triggered or sustained due to atrial drivers located elsewhere in the atria. We designed a new simulation workflow based on personalized computer simulations to characterize AF complexity of patients undergoing PV ablation, validated with non-invasive electrocardiographic imaging and evaluated at one year after ablation. We included 30 patients using atrial anatomies segmented from MRI and simulated an automata model for the electrical modelling, consisting of three states (resting, excited and refractory). In total, 100 different scenarios were simulated per anatomy varying rotor number and location. The 3 states were calibrated with Koivumaki action potential, entropy maps were obtained from the electrograms and compared with ECGi for each patient to analyze PV isolation outcome. The completion of the workflow indicated that successful AF ablation occurred in patients with rotors mainly located at the PV antrum, while unsuccessful procedures presented greater number of driving sites outside the PV area. The number of rotors attached to the PV was significantly higher in patients with favorable long-term ablation outcome (1-year freedom from AF: 1.61 ± 0.21 vs. AF recurrence: 1.40 ± 0.20; p-value = 0.018). The presented workflow could improve patient stratification for PV ablation by screening the complexity of the atria.


Introduction
Atrial Fibrillation (AF) is the most common arrhythmia, affecting a total population of 33.5 million worldwide [1]. Circumferential pulmonary vein isolation (CPVI) is considered the standard therapy for symptomatic AF patients [2]. However, non-pulmonary vein drivers located at the posterior wall, superior vena cava, the interatrial septum sites, the terminal crest or the coronary sinus can be found and are responsible in part for the inefficiency of the ablation procedure, especially in persistent AF patients [3].
Ablation planning and evaluation of atrial pro-arrhythmic behavior may play a key role toward ablation outcome. For this purpose, the combination of personalized atrial models implemented with electrophysiological and anatomical biomarkers of abnormal behavior have been integrated in a simulation environment to help identify arrhythmic behavior and improve novel diagnostic [4] and treatment strategies [5][6][7][8][9]. Computational simulations have emerged in this field as a new tool that can be used for characterization and prediction in different scenarios, from prediction of cardiotoxic compounds [10,11], targeted ablation [12] or recurrence after ablation [13].
In this field, automata models have been used for electrophysiological simulations to achieve simpler approaches with lowered computational time as compared to other models that include ionic level description. A lowered computational cost is translated into faster simulations with a higher number of possibilities to explore [14].
Here, we present a novel methodology to predict the efficacy of AF ablation based on computer simulations that included patient anatomy and different arrhythmic scenarios (i.e., different rotor location and number). These simulations were later compared with the clinical results of patients undergoing electrocardiographic imaging (ECGi) maps, CPVI and 1-year ablation outcome.

Patient Database
We included patients undergoing CPVI for drug-refractory paroxysmal (N = 14, 9 female) and persistent AF (N = 16, 8 female). Candidates were patients >18-years old, history of symptomatic AF, included if sustained AF was inducible during the electrophysiological study. Patients included in this study were admitted for ablation of drug-refractory paroxysmal and persistent AF, undergoing circumferential point-by-point ablation [2]. All patients gave informed consent. The protocol was approved by the institutional review board of the Hospital General Universitario Gregorio Marañón.

Atrial Electroanatomical Complexity Evaluation Protocol
Atrial electroanatomical complexity was evaluated analyzing the number and distribution of AF reentrant sites in relation to the anatomic characterization of the atrium. To that purpose: (1) MRI imaging from patients were obtained; (2) computational simulations of cardiac activity in the reconstructed atrium were performed; and (3) the results of simulations were compared with patients' clinical characteristics, ECGi complexity and outcomes after ablation.
The workflow followed for the evaluation of atrial pro-arrhythmic behavior is explained below and summarized in Figure 1.

Atrial Anatomy
Magnetic Resonance Imaging (MRI) was performed in all patients before ablation procedure. MRI images with a spatial resolution of 0.7 mm × 0.7 mm × 1.5 mm were acquired 2-3 days prior to the ablation procedure and segmented using ITK-SNAP [15]. Images were segmented to obtain a 3D mesh of both atrial cavities using growing region automatic segmentation for both atria separately. Later, Meshmixer software was used to combine both left and right atrium. After obtaining the raw anatomies, meshes were resampled to obtain a 200 µm resolution for simulations. An example of the final mesh used for simulations, which includes the atrial anatomical complexity present in patients, can be observed in Figure 1. Left atrial area (mm 2 ) was measured in every anatomy for further analysis together with the electrophysiological variables obtained from the workflow.

Computational Models of the Atria
Once the atrial anatomies were segmented, simulations were run under AF conditions where rotational activity could be characterized. Overall, for each anatomy, 100 simulations per patient run with different initiation protocols using the corresponding individualized anatomical model for 1000 ms. These simulations had an arbitrary location of rotational activity patterns on the atrial cavity ranging from 1 to 10 rotors simultaneously, ensuring the total coverage of the atria for later evaluation.
This protocol was repeated 10 times per geometry to increase the number of simulations, achieving a final number of 100 simulations per anatomical model. First, rotational activity was distributed over both cavities of the atria with different locations each time. After the location of the rotors was obtained [16], the automata model was run to evaluate the evolution of the scenario and posterior characterization. These models, despite having simpler formulations, allow for the presence of more complex scenarios, including location of higher number of rotors. The models for rotor location and activation model are explained in subsequent sections. A brief description of the difference between ioniclevel electrophysiological models and automata models showing examples in 2D planes is further discussed in Supplementary Material ( Figure S1).

AF Initiation: Automatic Rotor Location
Jacquemet et al. algorithm [16] was implemented for the development of automatic location of the rotational activity. This implementation, based on an eikonal-diffusion solver, allows obtaining computed activation maps similar to those obtained in the monodomain model, with the option of varying the number of rotors in the model from 1 single rotational foci up to 14 [16]. The location of the rotational activity was arbitrary and only depended on the curvature of the model. As Jacquemet's algorithm is calibrated in phase, a conversion into the labels for the automata model (3 discrete states) was performed to continue with the workflow. An example of this implementation can be found in Figure S2 together with the characterization of the simulations ( Figure S3).

Automata Model Simulations
An automata model based on activation patterns was implemented to simulate cardiac activity in the atrial cavity. Alonso-Atienza's model [14] was implemented to perform simulations with three different states (state 0 or resting, state 1 or activated, and state 2 or refractory period) that depended on the probability equation depicted below, where E represents the excitability of the unit, A the activation and D the distance matrix.
A more detailed description of this model, including all the equations and variables, can be found in the Supplementary Material. All nodes in the mesh were simulated following this model, i.e., no differences were implemented for different regions nor fiber orientation. NVIDIA Titan XP was used for all the simulations and posterior analysis of the workflow. Simulations were run in Microsoft Visual Studio 2017 and characterization of the simulations was performed in Matlab. The estimated ionic simulated model cost was 275 min vs. automata model: 42 min for 1-second simulation during AF, including stabilization and arrhythmia induction for the ionic model.

Electrophysiological Equivalence and Characterization
The evaluation of the electrophysiological properties of the simulations, which included the 3 states of the simulations of the automata, were calibrated using Koviumaki Action Potential Duration [17] to translate the automata model into measurable atrial electrophysiological signals. For this purpose, the square pulses that are identified as activations in the automata model, were directly substituted with the atrial APD morphology.
Once the electrophysiological information was recovered, electrograms were calculated for each node. More specifically, from each simulation, a uniform mesh of pseudounipolar electrograms was calculated under the assumption of a homogeneous, unbounded, and quasi-static medium [18]. The mesh used for the electrogram calculation was individualized and corresponded to the same mesh used for the ECGi calculation, allowing a direct comparison between both analyses.
In addition, the logarithmic energy entropy, which has been extensively used for the characterization of signals in other disciplines [19], as well as for cardiac signals [20], was calculated on the electrograms for each node and normalized for each atrial anatomy. More specifically, this entropy showed similar performance in prediction algorithms in previous studies [20] as Shannon entropy, widely used in the electrophysiological field. Finally, the mean entropy of the electrograms from all the simulations for a given patient was calculated and evaluated using entropy maps.
The main output of the workflow was produced by means of Atrial Complexity Maps (ACM) and Atrial Complexity Biomarker (ACB). ACM were obtained from the average entropy values of all the simulations from a given patient. ACB was obtained from the quantification of the number of rotors attached to the PV in the sustained simulations for each patient, which were later averaged. A rotor was considered to be attached if rotational activity was maintained around the PV for the complete simulation.

Clinical Evaluation AF Complexity: Atrial Complexity Map vs. ECGi
We compared the number of AF simulations with maintained reentries (ACM) obtained from the simulation workflow with the histogram of rotors obtained from the ECGi calculation. As explained in previous sections, the entropy maps were calculated with the same anatomies that the ECGi for them to be comparable. The specific protocol for obtaining and calculating ECGi was previously described [4,21,22]. Briefly, a minimum of three segments of at least 1 s duration were selected to calculate the histogram of rotors from ECGi signals. Rotors were obtained by counting the number of rotors in each atrial model node from the ECGi calculations for each of the segments. Finally, all three histograms were averaged and compared with the results of the simulations. Comparison was performed by dividing the complete anatomy by areas as shown in previous studies [23] and in Figure 2, and evaluating the presence of rotors and high entropy foci per area. This methodology enabled the characterization and evaluation of the complexity of the atria.

AF Complexity and 1-Year Ablation Outcome
Evaluation of the simulations was performed attending to the number of sustained simulated AF episodes per patient, rotors attachment to the pulmonary vein (ACB) and left atrial appendage, rotor distribution between both cavities, mean conduction velocity, and maintenance of more complex scenarios. These parameters were evaluated as predictors of ablation efficacy at 1 year. Additionally, AF type (Paroxysmal vs. Persistent) was also compared to reveal possible characterization patterns using the workflow.

Statistical Analysis
The t-test was used to evaluate the statistical significance between continuous paired or unpaired variables, and statistical significance was considered for p < 0.05 for continuous variables. Pearson Chi Independence test was used to evaluate categorical or binary variables, and statistical significance was considered for p < 0.05. All data are reported as mean ± SD. In addition, a regression analysis to test the independent predictive value has been conducted including the following parameters: AF type, gender, simulations results and 1-year outcome.

Cohort Description
In total, 30 patients were included in this retrospective study. The characteristics of patients included in the study are shown in Table 1. When patients were compared according to 1-year post-ablation outcome, there were no significant differences in age, height or weight. In contrast, the proportion of persistent AF patients at 1 year was significantly higher on the AF group and the same trend was observed for female patients.

Comparison of ACM with ECGI
In Figure 3 we show the coincidence between the histograms of rotors recorded in patients ECGi and high entropy areas obtained from the simulation protocol. Atrial Complexity Maps were identified as descriptors of the atrial complexity, where the characteristics observed on the histogram of rotors showed a direct correlation with 1-year post ablation outcome. This correlation showed 93.33% similarity in the pulmonary vein area and posterior wall, 80% coincidence in the floor area, 86.67% in the lateral wall and 83.33% in the right atrial appendage. An example of a simple Atrial Complexity Map is shown in Figure 4a where the entropy foci is mainly located on the left superior PV occurring in a patient that maintained SR at 1-year after ablation. Figure 4b shows a heterogeneous and complex Atrial Complexity Map with multiple high entropy foci, i.e., for electrograms that presented entropy values higher than 0.8*maximum entropy, distributed on both atria, with low rotor attachment to the PV area, which occurred in a patient with AF recurrence during follow-up.   The percentage of patients that presented high entropy values on the pulmonary vein area on the ACM was significantly higher on the freedom of AF group (n = 18) than of the AF group (n = 12) (freedom of AF: 93.75% vs. AF: 62.5%; p < 0.001), supporting the favorable ablation outcome in the former group. Moreover, freedom of AF patients presented lower number of high entropy areas or simpler ACM on the RA than AF patients (freedom of AF: 68.75% vs. AF: 100%; p < 0.001).

AF Complexity and 1-Year Ablation Outcome
Patients with freedom of AF at 1-year presented higher values of the Atrial Complexity Biomarker (i.e., a higher number of rotors attached to the PV) than patients with AF recurrence during follow-up (freedom of AF: 1.61 ± 0.21 vs. AF: 1.40 ± 0.20; p-value = 0.018) (Figure 4c). Interestingly, the mean number of sustained left atrial appendage rotors tended to be higher on the group of AF Freedom (freedom of AF: 3.52 ± 3.81%; AF: 1.50 ± 2.37%; p-value: 0.14) (Figure 5e). From the complete set of simulations, the number of rotors was higher on the left atrium (LA) than in the right atrium (RA) for both groups: AF Freedom patients (LA: 2.96 ± 0.35; RA: 2.22 ± 0.34; p-value < 0.0001) and AF patients (LA: 2.97 ±0.68; RA: 1.97 ±0.40; p-value < 0.0001) (Figure 5d). Differences in mean conduction velocity during AF were not significantly different between groups (freedom of AF: 85.91 ± 6.18 vs. AF: 85.73 ± 3.81 cm/s; p-value = 0.94) (Figure 5c). From all the possible simulated scenarios, some presented more stability in time, that is, a higher percentage of simulations sustained AF during 1000 ms. The relationship between the number of initiated rotors with respect to the number of rotors maintained in the simulations is shown in Supplemental Figure S3. In this case, the average number of sustained AF simulations presented a decreasing trend for an increasing number of initiated SP, that is, the arrhythmia was not easily sustained for high number of simulated rotors. Furthermore, there were no significant differences between the freedom of AF and AF group, except for the case in which the number of initiated rotors was fixed to two, in which the average number of maintained simulations of AF patients was higher in the AF group (freedom of AF: 3.87 ± 2.17 simulations vs. AF: 5.89 ± 2.09 simulations; p-value = 0.035). Further comparison of the electrophysiological description of the simulations can be found in Table 2.
Regression analysis results show that the ACM biomarker is the only variable with a trend for independently predicting 1-year post ablation outcome (p-value = 0.0752) as compared to other clinical variables such as AF type (p-value: 0.2548) and gender (0.3442).

Applicability to Clinical Environment
This methodology, and specifically, the ACM and ACB obtained from the workflow, are presented as an estimator of the complexity of the atrial activity. Figure 6 shows the rotor maps from three different patients ordered for increasing Atrial Complexity Biomarker values. As it can be observed, the higher the ACB, the lower the complexity of the rotor map with more localized the rotors appear on the pulmonary vein area. From left to right, patients with lower ACB presented higher number of rotors than patients with higher Atrial Complexity Biomarker. In addition, rotors in patients with low ACB were mainly located in both atria, whereas patients with high ACB presented the rotors concentrated on the PV area.
Examples of these cases are present in Figure 6, where the left panel, corresponding to an ACB with 1.39, represents a patient in which the ablation was not successful and both the ECGi rotor histogram and the Atrial Complexity Map present high activity foci on the right atrium. In contrast, the patient on the right panel represents a case in which the ablation was successful and both the ECGi rotor histogram and the Atrial Complexity Map present high activity foci on the right pulmonary veins, accompanied by a higher ACB. Overall, 11 patients (AF Freedom: 2 patients vs. AF: 9 patients) corresponded to an ACB lower than 1.45 and 19 patients (AF Freedom: 16 patients vs. AF: 3 patients) corresponded to an ACB higher than 1.50.

Discussion
This study presents a new simulation workflow for the personalized electrophysiological evaluation of the atria in a simulation environment. This is a proof-of-concept study that establishes a noninvasive evaluation of atrial electrophysiological complexity by means of two novel biomarkers: Atrial Complexity Maps and the Atrial Complexity Biomarker to evaluate atrial complexity and predict the efficacy of the ablation. Our results revealed that long-term successful AF ablation occurred in patients with rotors mainly located in the pulmonary veins, while unsuccessful procedures presented greater number of entropy foci outside the PV areas. Paroxysmal AF patients presented a significantly higher number of LA rotors with greater attachment to the PV area and lower density of entropy areas in the RA, as compared to persistent AF patients, explaining the improved long-term clinical outcome in the former group. These results are in agreement with clinical studies that consistently report a better outcome of the ablation following PV electrical isolation, while AF recurrences were most common in patients with multiple drivers distributed at extra-PV sites. [2,3,21,24,25].

Simulation Models
Detailed ionic models for the evaluation of AF induction at different scales include complementary information that may be relevant for the ablation procedure [11,26]. These ionic models also appear as a good approach for very specific workflows, such as pharmacological studies, which analyze the playing role of ionic channels and its modification using different drugs. In addition, several studies explore different strategies for the evalu-ation of the atrial complexity, including specific ablation strategies and targets [5], cycle length evaluation [27] and proarrhythmic structures such as the left atrial appendage [28]. However, they may present significant drawbacks, especially for large-scale simulations, due to the high computational power needed for such specific models, limiting the overall number of scenarios to be studied [29], therefore limiting the number of scenarios to be studied or the number of structures (i.e., only including left atrium). In contrast, simpler models, such as the automata model used in this study, can be applied for modeling an initiated arrhythmia behavior enabling the analyses of several distributions of rotational foci. Moreover, automata models rely on simpler activation patterns, and can be implemented and used on cardiac modeling to obtain similar approaches with a lower computational cost [14,30]. Therefore, the use of simpler models together with graphical processing units for parallel computation, reduces the total computational time, allowing a potential translation and implementation of this methodology in the clinical environment for patient evaluation.
These simulations are presented, as a workbench for characterizing the proarrhythmicity based on the anatomy and different arrhythmic scenarios. One of the main challenges in computation is the initiation of rotational activity on the desired area. Several approaches have been implemented and described in previous publications in order to tailor arrhythmia initiation by including remodeling such as repolarization alternants, adipose tissue modeling, and cardiac ion channel mutations [8]. However, we gave priority to the analysis of scenarios with different combinations of rotational activity that reflect the heterogeneity of the arrhythmias using an algorithm that directly deal with different rotors over the atria, comparing their different distributions. The inclusion of such a high number of scenarios or combinations of rotational foci (i.e., 100 simulations per anatomy) enables to include all possible areas at which rotors can be maintained, differing from other approaches in which a low number of combinations is analyzed, restricting the arrhythmic simulations to the pulmonary vein area and excluding the arrhythmia initiation on right atrium [31].
Regarding the characterization of the simulations, all simulated atria presented realistic models in which the number of rotors was higher on the left atrium than in the right atrium, with a similar number of maintained simulations per group and high attachment of rotational drivers to the pulmonary vein area, identified as the main proarrhythmic trigger on clinical practice. These results align with previous studies that reflect the dominance of the LA in the rotational activity of AF patients [4,24,25,[32][33][34][35], demonstrating the reproduction of a clinical scenario into personalized simulations in a computer.

Clinical Implications
The increasing number of potential candidates for ablation therapies is much higher than the availability of laboratories to perform procedures, but patients are selected based on very simple and unproved selection criteria for efficacy. However, current indiscriminate application of ablative therapies to large, unselected cohorts of patients with atrial fibrillation might dilute the intended treatment benefits and significantly increase the cost. Translation of the mechanistic insights of computational and basic research into clinical management concepts will uncover the full potential of personalized atrial fibrillation management [36]. The present approach may help appropriately select patients undergoing invasive therapies by: (1) integrating the workup protocol as shown in Figure 7, where anatomical characterization and simulations will be performed 2 days prior to the procedure, to later correlate high entropy areas location in the simulations protocol with ECGi, and help decide the ablation strategy; (2) giving preference for the standard ablation procedures to those patients with favorable predictors for ablation long-term success (low ACM, high ACB) [37]; and (3) selecting patients with higher atrial complexity to undergo the elimination of extrapulmonary drivers ablation [21,32,38]. This protocol is based on a personalized simulation method that could be potentially modified to input other remodeling factors such as fiber orientation to develop more complex fibrotic models, broadening the application to models closer to the clinical setting. Therefore, the integration of image-based computational modeling into treatments for heart rhythm disorders could thus advance personalized approaches to heart disease.

Study Limitations
The main advantage of the model-its simplicity-also constitutes its main limitation, as this model is not as tailored as ionic models. Furthermore, other scenarios such as different conduction velocity areas or fiber orientation should be implemented on the automata model and considered to study a wider population of patients. Second, we were able to demonstrate that all the areas that presented high-frequency activation on the ECGi presented high entropy values, but we were unable to ensure that all the high entropy areas were coregistered with high-frequency areas or rotational activity. Further studies should be conducted to evaluate if this mismatch is due to a lack of more ECGi episodes or if rotors were present only in part of the atrial anatomy.
In addition, other tailored characteristics such as fibrosis distribution over the atria and ionic heterogeneity [39] should be considered in further studies to better represent the anatomical and electrical remodeling of the cardiac tissue. Atrial thickness and blood pressure are two important factors that have been demonstrated to affect frequency dynamics and should be further explored to complement the models [40].
Finally, higher attachment to the left atrial appendage was observed on the AF freedom group, which exclusively underwent PV isolation. Further studies should confirm the proarrhythmic behavior of the left atrial appendage in these models [41].

Conclusions
This study presents a new method for the evaluation of the pro-arrhythmic areas on atrial anatomies providing Atrial Complexity Maps and the Atrial Complexity Biomarker as estimators of atrial complexity. This approach, validated using ECGi to measure atrial complexity, was able to identify the set of patients that presented higher atrial complexity. Informed Consent Statement: Due to the retrospective nature of the study, informed consent was waived by the IRB.