- 1Department of Physics and Astronomy, Ghent University, Ghent, Belgium
- 2Department of Pharmaceutical Sciences, Antwerp University, Antwerp, Belgium
Introduction: Atrial fibrillation (AF) is a very common cardiac arrhythmia whose mechanisms are still a topic of debate. This work applied concepts of topology to gain new insights into reentry-based simulated AF, similar to our previous work in atrial tachycardia (AT). We demonstrate that the Index Theorem – which states reentries must come in pairs of opposite rotations – applies to a model of AF, even when the complex dynamics change over time. Additionally, we tested the hypothesis that connecting opposite pairs of singularities can terminate simulated AF in the same way as clinical and simulated AT.
Methods: we applied a modified phase mapping capable of detecting both functional and anatomical reentry to a dataset of 600 AF simulations based on clinical data. We then compared three virtual ablation strategies: random lines, straight lines, and heuristic lines. Straight lines connected pairs of opposite singularities through the shortest path; heuristic lines connected them in such a way that prioritized blocking the conduction path; and random lines connected randomly selected pairs of points with comparable distance to the other methods.
Results: We showed that our algorithm could verify the predicted paired reentries for 99% of the simulation duration on average, and 93% for the worst-performing case. The heuristic virtual ablation method terminated activity for 90% of cases, a marked improvement against the straight line method (55%) and the random method (0.5%).
Discussion: This work provides mechanistic insights into AF, and points towards pitfalls of ablation strategies, both of which have the potential to improve our understanding and ability to treat this condition.
1 Introduction
The mechanisms behind atrial fibrillation (AF), the most common cardiac arrhythmia, are still not fully understood. Here, we apply concepts from topology to reentry-based simulated AF, expanding our recent works in which a similar methodology brought new mechanistic insights into atrial tachycardia (AT) (Duytschaever et al., 2024; Abeele et al., 2025).
In the context of cardiac electrophysiology, phase mapping associates the different states a cardiomyocite or patch of cardiac tissue may assume to different values in a limited range. Typically that range is
For a given phase field
A phase singularity, or singular point, is defined as a point around which any sufficiently small curve
In our previous works, we extended the concept of phase singularities, showing that the Index Theorem holds not only for singular points, but also when integrating around anatomical boundaries. As a consequence, anatomical reentries bounded to a two-dimensional surface must always come in pairs of opposite chirality, represented by an index of
This fundamental property had been previously overlooked, largely due to the presence of incomplete reentries–circuits that fail to complete a full rotation before colliding with another wave. Historically, these passive reentries were considered unimportant for the maintenance of AT.
However, we demonstrated that an incomplete reentry–sometimes called a passive, bystander, or wannabe reentry (Lee et al., 2015; Maury et al., 2019) – is just as critical for sustaining AT as a complete (or full, driving, active) reentry. Although it does not complete a full rotation in physical space, the line integral around it still results in a non-zero value, describing a full rotation through phase space. Consequently, if ablation only targets the more easily identifiable full circuits while neglecting incomplete ones, arrhythmia will not terminate, and the incomplete circuit may even evolve into a complete one. Therefore, both complete and incomplete reentries must be ablated to achieve termination. Due to the Index Theorem, each clockwise rotation is matched by a counterclockwise counterpart, which may be of either type.
Based on this understanding, we classified anatomical boundaries between 3 types: critical boundary type 1 (CB1), which hosts a complete reentry; critical boundary type 2 (CB2), which hosts an incomplete reentry; and non-critical boundary (NCB), which does not host a reentry. Despite clinical terminology often referring to complete reentries as “active” or “driving” reentries (Santucci et al., 2024), any combination of CB1 and CB2 pairs may be present in a given case, depending on geometry and distribution of conduction velocity values (e.g., CB1-CB1, CB1-CB2, CB2-CB2). Regardless, all CBs must be connected by an ablation line to a CB of opposite index in order to terminate the arrhythmia. All three categories are illustrated in Figure 1. Note that multiple waves may interact with the same boundary, and the net total in each direction will determine the corresponding index, as described previously.
Figure 1. Diagram representing the previously described interactions between propagation waves (arrows) and a boundary (in gray). Each arrow’s head marks a wavefront, and its trailing line marks a waveback. From left to right, CB1 (complete reentry), CB2 (incomplete reentry), and NCB (no reentry). Each wavefront contributes to the index by a factor of +1 or −1 depending on direction.
Now, we seek to verify the Index Theorem for simulated AF, and explore its mechanistic implications. A commonly proposed AF mechanism is the presence of not only anatomical reentry but also functional reentry (rotors). According to the Index Theorem, singularities must come in pairs, including both functional and anatomical reentries. Paired rotors of opposite chirality have been described in multiple studies (Xu et al., 2023; Atienza et al., 2015; Gurevich and Roman, 2019). However, paired rotors are often described as a mere alternative to single rotors, which are considered an equally valid configuration (Lin et al., 2013; Narayan et al., 2014; Nattel et al., 2017; Rappel et al., 2024). That is likely because those works did not consider other types of singularity these rotors may be paired with.
In theoretical works focused specifically on singularity topology, it has been demonstrated that the topology of electrical propagation requires that phase singularities are only created or destroyed in counter-rotating pairs (Marcotte and Grigoriev, 2017; DeTal et al., 2022; Arno et al., 2024). They show that the total sum of the index, or topological charge, is a conserved quantity, analogous to the Index Theorem. However, boundary interactions were either not analyzed, or treated as a problem that can break this conservation. This limitation weakens the value of this topological approach, as wave collisions with boundaries are common.
Here, we reconcile our previous work, where rotation around boundaries is shown to obey the Index Theorem, with the above studies, which overlooked the contribution of boundaries when calculating the total topological charge. Topologically, an anatomical reentry is equivalent to a non-meandering rotor, as the boundary can be continuously reduced to an arbitrarily small core, like that of a rotor. Alternatively, the refractory core around which a rotor is anchored, or around which it meanders, can be treated identically to an anatomical reentry (Arno et al., 2024; Tomii et al., 2021). Moreover, rotor tips are topologically equivalent to other phase singularities with no special name, associated to the endpoints of wavefronts (Marcotte and Grigoriev, 2017) and non-excitable cores (Arno et al., 2024; Tomii et al., 2021). Any singularity has the potential to become a rotor tip, depending mainly on its interactions with other waves and refractory tissue, and whether those interactions result in sustained reentry or other behavior, such as meandering.
Therefore, as we have already showed in previous works that the Index Theorem holds for regular, anatomical reentries, and there is a topological equivalence between anatomical and functional reentries, it is to be expected that it also will hold for irregular, functional reentries–which we investigate in this work.
This work aims to show that the Index Theorem holds for simulated AF data that includes a mix of functional and anatomical reentries, often with unstable activity, waves breaking and merging, and boundary interactions. Furthermore, we show that the identification of singularities is fundamental for AF termination through virtual ablation. As the previous related works (Marcotte and Grigoriev, 2017; DeTal et al., 2022; Arno et al., 2024) were more concerned with the theoretical aspects, they only show a small number of examples, which were not anatomically realistic and did not consider boundary effects. To demonstrate it at a larger scale, we performed 600 simulations based on a previously described set of 100 anatomically modeled 3D meshes, which include realistic fiber direction, MRI-based fibrosis levels, and remodeling. Moreover, expanding on our previous work, we compare different ablation methods, obtaining unique mechanistic insights into the dynamics of irregular activity.
2 Materials and methods
To demonstrate our concepts, we generated a set of 600 AF simulations, from which we extracted the phase from the action potential signals. For each time-step, we identified and clustered phase singularities. Additionally, for select points in time, we used a modified distance based on the phase value to consider possible ablation lines between opposite-chirality cluster pairs; and finally, applied the combination of ablation lines that minimized the total ablation length. As a baseline for comparing outcomes, for each simulation we also generated a set of geodesic lines connecting opposite-chirality pairs, and a set of random ablation lines with approximately the same length.
2.1 Simulated data
Single Cell Model: All simulations were executed using the OpenCARP simulation software (Plank et al., 2021). We used the Courtemanche ionic model with AF remodeling as the electrophysiological model of our simulations (Courtemanche et al., 1998). This remodeling encompasses a reduction of the ultra-rapid delayed rectifier conductivity
Rather than modeling fibrotic regions as non-conductive tissue interspersed with conductive tissue, they were treated as conductive tissues with a further adapted ionic model.
These parameter changes resulted in a 14% reduction of CV relative to non-fibrotic regions.
Transmembrane voltages were recorded in 5 ms time-steps (200 Hz sampling frequency).
Mesh parameters: For our simulations, we used a publicly available dataset by Roney et al. (2022) consisting of 100 3D meshes of left atria based on clinical AF data, including MRI intensity and fiber orientations. It included 43 paroxysmal, 41 persistent, and 16 long-standing persistent AF patients, none of which had previously gone through ablation treatment. The meshes were refined to achieve an average edge length of 160–170 μm for better simulation precision. To mimic fibrosis, different conduction velocity values were assigned based on MRI intensity, as described in Table 1. Conduction velocity was then further reduced by 50% in the directions orthogonal to the fiber orientation.
Reentry induction: On each mesh, we induced reentry by means of phase singularity distribution. We distributed phases from -
Figure 2. Example of the 6 starting conditions for the same mesh. Singularities in opposite directions are placed along each of the Cartesian axes (from left to right: x, y, and z axis). The top row shows the clockwise placement, and the bottom row shows the counter-clockwise placement.
All simulations were analyzed starting from time = 50 ms until the final time of 3,400 ms, or until early termination, whichever happened first. A simulation was considered to terminate early, i.e., to spontaneously terminate, if all activity ceased without the need of any ablation. The termination time was defined as the moment after which all transmembrane voltages remain below −40 mV.
2.2 Phase singularity detection
While a fine-resolution 3D mesh was necessary to obtain accurate dynamic behavior for simulations, our subsequent analysis could be performed with coarser resolution, reducing computation times significantly. Therefore, apart from simulations, all analysis was performed on a subsampled version of the mesh in order to obtain approximately even-spaced points at a distance of 2 mm from each other.
To identify phase singularities, we applied a modified phase mapping method. We obtained local activation times (LATs) from the transmembrane voltages of the sampled points, and calculated their phase
Figure 3. Example of an original simulation mesh (left) and reduced mesh (right). The reduced mesh displays the geodesic Voronoi diagram connecting neighboring sampled points with white edges, and includes an additional point at the center of each anatomical hole, allowing the identification of anatomical reentries around it.
As explained in the Introduction, counting the number of discontinuities caused by the presence of wavefronts is a simple way to obtain the phase index. To do so, each point’s neighbors were first sorted in counterclockwise order, and phase differences between consecutive neighbors were computed. A phase difference was classified as a phase jump if
Points with an unequal number of phase jumps in each direction had a nonzero index and were tagged as singularities at that time step. This definition accommodates multiple jumps in each direction. In the terms of our previous AT works (Duytschaever et al., 2024; Abeele et al., 2025), a complete circuit (CB1) will have just one phase jump for its entire duration, while an incomplete circuit (CB2) will have two jumps occurring in one direction and one jump in the other for some of its duration. It is also possible for a boundary to have index values
After being identified, singularities were clustered, first in space and then in time. At each timestep, given the mean distance between neighbors in the subsampled mesh
This approach allowed us to identify which detections belonged to the same singularity and to measure each cluster’s lifespan, accounting for meandering. In all metrics described below, each cluster was treated as a single singularity rather than multiple detections. For example, a cluster of three points with index of −1 each, only contributed −1 to the total index, not −3.
2.2.1 Adherence to the Index Theorem
We defined two performance metrics based on our algorithm’s ability to detect singularities and track them over time–and therefore they are measures of the algorithm’s adherence to the expected results for the Index Theorem.
They are two relative measures:
2.2.2 Arrhythmia complexity
To verify that our simulations were suitably complex for AF, we estimated how long singularities lasted, the total number of singularities over time, and the highest number of simultaneous singularities. Unlike the Index Theorem consistency measures, we did not make a distinction between positive and negative singularities here. Between these three measures, we can observe for each simulation the stability of its singularities, how often they tended to be created and destroyed, and the largest number of simultaneously hosted singularities.
For a more clinically grounded perspective on the data, the cycle length (CL) median and interquartile range were calculated for all simulations and compared to typical clinical values. The activation times of each vertex were calculated from the transmembrane voltages, and the differences between consecutive activations obtained. Then, using the distribution of activation time differences for a simulation, the median and interquartile range–the difference between the 75th and 25th percentiles–were calculated.
2.2.3 Comparison of starting conditions
Since each of the 100 meshes were simulated under 6 different starting conditions, it is relevant whether there were relationships between simulations originating from the same mesh. Therefore, we recorded how often the simulations from each mesh terminated early and whether singularity hot spots were spatially correlated between different starting conditions for the same mesh.
To check for singularity hot spots for a given simulation, we calculated the singularity count
To compare the distribution of
2.3 Ablation strategy
To obtain mechanistic insights on the effects of ablation in our simulations, we generated virtual ablations. These are simulated lines of conduction block created instantaneously, idealized approximations of a real ablation. Different strategies were used to determine the size, shape, and position of each virtual ablation.
Heuristic lines: Virtual ablation lines were created for each simulation at the points in time of 1,000 ms, 2000 ms, and 3,000 ms. At each time-step, points within each cluster were treated as one singularity, and if a cluster contained multiple points, a conduction block loop was created passing through all of them. Virtual ablation lines were then created to connect pairs of opposite singularities while blocking the wave path and minimizing virtual ablation length. If there was a detection issue that caused an unequal number of opposite singularities, we iteratively increased the detection time-step by 5 ms until the index sum was 0.
To create a virtual ablation line between a pair of opposite-chirality singularities, the Dijkstra traversal algorithm was used to find the lowest-cost path between the two points. In order to obtain the desired effect of blocking wave paths while minimizing distance, a heuristic distance
As
The path calculation was repeated for each possible combination of opposite-chirality singularities, and the set of combinations that resulted in the smallest total virtual ablation (in heuristic distance) was chosen and performed. Because of the heuristic, these virtual ablations tend to be between singularities that belong to the same wavefront, but on occasion may connect points from different waves–for instance, if the distance between them is small enough to overcome the penalty to not ablating at the wavefront. Therefore, the method prioritizes blocking waves directly, but allows forcing collisions between different waves when singularities are sufficiently near each other.
Straight lines and random lines: As baselines to compare the heuristic virtual ablation lines, we generated two other types of virtual ablation lines for each simulation at the designated virtual ablation times: straight lines and random lines. In the straight line method, we connect singularities using the same logic described above, but using the geodesic distance over the surface rather than the heuristic distance–the shortest geometric path in the manifold of the reduced mesh. In the random method, we generated straight lines of comparable size to our heuristic virtual ablations, but connecting random pairs of points rather than singularities. To do so, we calculated the geodesic distance between all pairs of points in the subsampled mesh, and divided the pairs into distance deciles (10-percentile bins). For each singularity-based virtual ablation line, a new pair was randomly selected from the same distance decile as the corresponding original pair, assuring the random lines have comparable size. Figure 4 illustrates the three types of virtual ablation lines.
Figure 4. Example of the three types of virtual ablation line, from left to right: heuristic line (pink), straight line (blue), and random line (orange). Both the heuristic line and straight line form a loop around the same singularity detections and then connect them, but the former blocks the path of the wave, while the latter takes the shortest path. The random line connects the shortest path between two random points at a comparable distance to the other virtual ablation lines.
Once a virtual ablation line was created using the subsampled mesh, it was translated back into the full mesh by finding geodesic paths between consecutive sampled points. This way, the state of the simulation at the moment of virtual ablation could be re-created, and block lines placed instantly at the corresponding point in time. A virtual ablation was deemed successful if all activity ceased 400 ms after creation of the block line. For the virtual ablation simulations, transmembrane voltages were recorded in 10 ms time-steps (100 Hz sampling frequency).
2.3.1 Spontaneous termination
Virtual ablations were still performed on simulations that terminated spontaneously, but no virtul ablation was performed if one of the virtual ablation times (1,000 ms, 2000 ms, 3,000 ms) was after the termination time, as there would be no effect. However, virtual ablation was still performed for times before the termination time–e.g., if the early termination time was 2,500 ms, the virtual ablations at 1,000 ms and 2000 ms would be performed, but not the one at 3,000 ms. We recorded which simulations terminated early, what their termination times were, and how virtual ablation affected termination time.
Additional criteria were required to determine virtual ablation success for spontaneous termination cases. Each outcome was simulated from the moment of virtual ablation
• If
• If
• If
3 Results
3.1 Representative examples
Before showing our general results for the dataset, we illustrate the TC detections and clustering with some representative examples. Figure 5 shows a snapshot of a case in our dataset where the index sum was 0 for its entirety, at a time when there were 4 positive and 4 negative singularities. The color-coding highlights how each wavefront is associated with a pair of opposite singularities, and that these pairs may be either functional-functional or functional-anatomical.
Figure 5. Snapshot of a successful detection, viewed from different angles. The mesh is colored according to phase values. Phase jump detections are marked with orange arrows in the direction of propagation. Detections of positive and negative singularities are marked by red and blue dots respectively. A label is placed next to each cluster of points, indicating that cluster’s contribution to the total index, and matching the color of the opposite-direction singularity associated with the same wavefront. Note that the pair marked in green is an anatomical reentry (mitral valve) paired with a functional reentry.
Figure 6 (top panel) illustrates the time evolution of the number of singularities and index sum for our worst-performing simulation in terms of the measure
Figure 6. Top: Plot of the number of positive and negative singularities (solid red line and dashed blue line, respectively), and the index sum (solid black line), over time for a given simulation. Positive and negative singularities vary over time, but the index sum remains 0 for the majority of the time, with small fluctuations. Bottom: Snapshot of a moment when the index sum did not equal 0, viewed from different angles and represented with the same conventions as Figure 5. To the right, the image zooms into a section of the simulation where an incorrect detection occurred. Circles highlight which singularity detections were assigned to which clusters in the zoomed section. An unpaired −1, marked in red, results in a non-zero total index.
To demonstrate the adherence to the Index Theorem, we used
The bottom panels of Figure 6 show a snapshot where the index sum did not equal 0 for that simulation. Note how most singularities are paired off as expected, with a wavefront associated to each pair, but an extra
3.2 Phase singularity detection
3.2.1 Adherence of the Index Theorem
Figure 7 shows histograms of our two performance measures across all simulations of the dataset. All simulations had a very high adherence to the expected Index Theorem result, with the lowest value of
Figure 7. Histograms for the two performance measures. Left:
3.2.2 Arrhythmia complexity
Figure 8 shows the distribution of values for our three complexity measures. On the left, we see that the majority of singularities are very short-lived, with a median value of 30 ms. However, some outliers are significantly more stable, persisting for most of the simulation–in one instance, a singularity lasted for all 3,400 ms of the simulation. On the center, we have the distribution of total singularity counts for the entire duration of each simulation. While the most stable cases had only 2 singularities for their entire duration, in the majority of cases many singularities are created and destroyed over time: a median count of 62 singularities over the simulation, and some cases reaching upwards of 250. On the right, the distribution of maximum simultaneous singularities shows again that the simplest cases had only 2 simultaneous singularities, but the median case had 10 and some had up to 26.
Figure 8. Histograms for the three complexity measures: Left: Singularity duration. While a few singularities last for most of the simulation, by far most of them are rather short-lived. Center: Total singularity count over time. As singularities are created and destroyed through the evolution of each simulation, many singularities can be detected. Right: Maximum simultaneous singularity count. Unlike AT, where only two counter rotating reentries are typically observed, here we see up to 26 singularities co-existing in some cases.
Figure 9 shows the distribution of values for our CL measures. Clinical AF CL values range from around 90 ms–250 ms and are typically centered around 180 ms (Haïssaguerre et al., 2004; Nagy et al., 2021). Our distribution of median CL, while within that margin, leans towards slower values. However, the spread associated with the interquartile range suggests some cases had significantly faster regions or moments of activity, consistent with the irregularity of AF.
Figure 9. Histograms for the three complexity measures: Left: Distribution of median CL values. Right: Distribution of CL interquartile range values.
3.2.3 Comparison of starting conditions
For each mesh, we calculated the Moran’s
Figure 10. Histograms for the two Moran’s
3.2.4 Spontaneous termination
29% of the dataset spontaneously terminated. Figure 11 shows on the left the distribution of early termination times across those cases, with a median of 1952 ms. Figure 11 shows on the right that depending on the mesh, anywhere between 0 and 5 out of 6 starting conditions resulted in early termination, though most meshes leaned towards not terminating early. Notably, all meshes had at least one starting condition that could sustain activity for the full duration rather than spontaneously terminating.
Figure 11. Histograms for the two termination measures. Left: Distribution of early termination times. A wide variety of termination times were observed, with some lasting almost the full duration of the simulation and some terminating before even the first virtual ablation time. Right: Distribution of early termination counts–how many starting conditions terminated early per mesh. 23 meshes never terminated early, and only 3 terminated early for 5 out of 6 starting conditions. Other meshes fell somewhere between, but notably no mesh terminated early for all 6 starting conditions.
3.3 Virtual ablation outcomes
As previously described, virtual ablations were not performed when a case spontaneously terminated before the chosen virtual ablation time, resulting in a total of 1,617 simulations per virtual ablation strategy (out of a possible 1800, had no early terminations occurred). Following the success metrics for each of the virtual ablation strategies, we separately calculated the performance for the spontaneous termination cases and for the non-spontaneous termination cases, as well as the overall performance, as seen in Table 2. The results are consistent between the two groups, except for a minor increase in the very low performance of random lines. The straight lines showed a marked improvement compared to the baseline of the random strategy, but still failed roughly half the time. The heuristic lines on the other hand, almost always succeeded, though they performed better for the non-spontaneous cases.
4 Discussion
4.1 Paired rotations and ablation in AF
In this work, we extended our previous results from AT to simulated AF, demonstrating how the Index Theorem holds for 600 patient-based AF simulations. Note that algorithmic failures to detect an index sum of 0 do not suggest that the Index Theorem is not true: they show that even with an imperfect algorithm, it is possible to observe its adherence for the majority of the time, reinforcing its relevance for AF dynamics. We therefore generalize our previous results for both functional and anatomical reentries, as well as irregular activity.
We also showed that opposite-chirality singularities must be connected to terminate the arrhythmia. However, we observed that only connecting the singularities with a block line is often insufficient, and that the particular way that these connections are made plays a major role. The necessity of connecting all singularities suggests that arrhythmic activity is sustained not only by long-lived anatomical reentries and rotors, but also by shorter-lived singularities, as ignoring them during virtual ablation may lead to continued activity.
As the heuristic method was designed to block the wave paths, its high performance is expected, and its significantly higher value compared to the straight line method highlights the significance of the particular way singularities are connected by ablation lines–in contrast to our previous AT work, where simply connecting pairs of CBs was sufficient. The heuristic method was a “hard limit”, in which the path was completely blocked; but less aggressive strategies that still block conduction paths and do not allow reoccurrence may be possible.
Observing the 8% of failed heuristic virtual ablations in cases that did not terminate early, we noted that all of these cases had something in common: despite being designed to do so, the conduction block lines did not fully block the waves at the moment of virtual ablation. This arises from the specific design of our heuristic, which weighs paths between points using a combination of their geometric distance and phase values at that point in time. Depending on the relative sizes and positions of different waves, these two weight components may act in opposition: large geometric distances increase the weight, whereas close phase alignment with a wave decreases it. In such cases, the minimum-weight path may connect singularities from different waves, thereby failing to block the waves as intended. The fact that this happened in every case where the heuristic lines failed strongly points towards the idea that those failures were caused by improper connection of singularities, again suggesting that very particular forms of ablation are required to truly terminate the arrhythmic activity.
As supporting results, we also analyzed the simulation dynamics, showing that the simulations were fittingly complex for AF, outside of a small number of outliers. This was a particularly relevant concern due to the relatively artificial starting conditions, based on phase singularity placement. However, as we observed, tens or even hundreds of singularities were created or destroyed per simulation. Most of these singularities were rather short lived, and most cases hosted multiple of them simultaneously, demonstrating the evolution away from the simplified initial setup. Almost 30% of simulations terminated spontaneously, but all meshes had at least one starting condition that did not. Notably, not only was inducibility greatly affected by the starting conditions, but also very little spatial correlation was found between singularity hot spots from different simulations on the same mesh–possibly suggesting a secondary role of the substrate, with starting conditions being more significant, as is expected of a chaotic system.
4.2 Paired reentries in the cardiac arrhythmia literature
While previous theoretical studies confirmed that the Index Theorem applies to AF (Marcotte and Grigoriev, 2017; Gurevich and Roman, 2019; DeTal et al., 2022; Arno et al., 2024), they had limitations that we address here, generalizing their findings. The authors in Marcotte et al. (Marcotte and Grigoriev, 2017) demonstrated that wavelets must be bounded by pairs of opposite singularities, but did not analyze how the singularities would interact with non-conductive boundaries. Gurevich et al. (Gurevich and Roman, 2019) claimed the total topological charge is conserved, except during boundary interactions–however, as showed here, by computing the boundary’s topological charge, this inconsistency is easily solved. By generating additional wavelets to induce collisions between singularities of opposite chirality, DeTal et al. (2022) were able to effectively eliminate all rotational activity. While still considering the existence of “single” rotors, they demonstrated that such rotors could be eliminated by treating the boundary as a mirror, effectively connecting each rotor to its “mirror image” with opposite chirality. This approach, while effective for generating opposite wavelets, still assumes that the total topological charge is not conserved when dealing with boundaries. Here, we show that, in fact, by simply considering obstacles as also being able to host reentries, the inconsistency is easily patched, and the Index Theorem is preserved. Arno et al. (2024), somewhat as a side result, showed that the Index Theorem must hold, but did not examine boundary interactions, as their interest was largely on meandering rotors with linear cores.
Moreover, because all of these works are focused on theoretical developments, their results, both simulated and experimental, are shown in rather small datasets, seeking to illustrate more than demonstrate. The data itself was either simulated in 2D tissue or experimentally observed from cell cultures. In this work, we sought to overcome this with a large dataset that incorporates clinically-derived anatomical and structural information, showing both functional and anatomical reentries can be treated in the same way.
Perhaps even more notably, it is still common in clinical work to describe the presence of both “single” and “paired” rotors (Lin et al., 2013; Narayan et al., 2014; Nattel et al., 2017; Rappel et al., 2024). Although theoretical knowledge has existed for years, albeit with limitations, it has not yet been widely spread in the clinical side of the field. We noted the same in our previous AT work (Duytschaever et al., 2024; Abeele et al., 2025), with clinicians often not making ablations that connected all singularities, and therefore unknowingly risking the possibility of arrhythmia reoccurrence. While our own simulated ablation strategy is not realistically applicable, it gives an indication towards the paths that should be taken in future developments.
In our previous AT work, we identified boundaries as of 3 types: NCB, CB1, and CB2 (Duytschaever et al., 2024; Abeele et al., 2025). We correctly identified that CB1 and CB2, the boundaries that host reentries, must be connected by an ablation line to terminate the reentrant arrhythmia. With the context of this work, we can refine that understanding. As we showed here, simply connecting opposite singularities does not guarantee termination, and is highly dependent on the particular choice of connection. While it eliminates the currently present singularities, the possibility that new singularities may emerge later means that the arrhythmia may reoccur. In our previous work, as long as all opposite singularities were connected, we could expect termination–thanks to the characteristic regularity of AT, which guarantees that eventually all waves would collide against the ablation lines and have nowhere to propagate, terminating all activity. This continues to be true, as demonstrated in initial clinical trials (Duytschaever et al., 2024), but it is a property of AT that cannot be extended to AF. That is why in this work we tested both the geodesic line virtual ablation, and the heuristic wavefront-based virtual ablation, showing that termination is only certain if conduction block not only connects two opposite singularities, but also makes it impossible for new singularities to arise later. Our chosen approach was a “brute force” upper limit to demonstrate the concept, but future work may investigate a way to optimize ablation configurations, and possibly test more realistic implementations.
4.3 Limitations
Our dataset consisted exclusively of simulated data, as our analysis would require global simultaneous AF data at a resolution that is not currently clinically available. Future works may use different models, or experimental data, to further verify our observations.
The presence of functional reentries in clinical practice remains controversial, as they have so far been primarily reported in computational studies (Xu et al., 2023; Allessie and de Groot, 2014; Waks and Josephson, 2014). Consequently, our findings are mainly relevant in the context of AF maintained by rotational activity. Nevertheless, because most computational studies on AF describe reentries as the dominant sustaining mechanism, our work directly contributes to this line of research (Narayan et al., 2014; Morgan et al., 2016; Vigmond et al., 2016; Haissaguerre et al., 2014).
As is common in simulation studies, we aimed to ensure physiological realism by incorporating MRI-derived atrial geometries and fiber orientations, and by employing the Courtemanche model, which captures key electrophysiological properties of human atrial tissue. Within this framework, we observed the emergence of multiple reentries, both stable and unstable.
Our chosen induction method, based on OpenCARP’s phase distribution approach, is less realistic than stimuli delivered by either the cardiac conduction system or from pacing during treatment. However, as we showed with our complexity measures, the majority of simulations evolved dynamically over time in organic ways that distinguish it from the initial setup.
Our virtual ablation strategy is not clinically realistic, as we instantly create conduction block lines, which in a real setting would move over the time needed to ablate. Moreover, our observations relate only to acute termination, and we have not verified re-inducibility or long-term effects. However, that was not the focus of this study. Rather, this is meant as a proof of concept to demonstrate the importance of the Index Theorem for our AF model, and our assertion of how connecting singularities is necessary, but not sufficient. It is not a viable recommendation, but rather a demonstration of issues with the current understanding of how to achieve termination through ablation.
As previously mentioned, while having a high performance, our algorithm could not detect an index sum of zero for all times. This was largely due to issues with clustering, and brief failures to detect phase jumps. While these errors do not significantly affect our observations and conclusions, we aim to introduce a future work with an updated version of the method that identifies singularities more robustly.
Our simulations cover common proposed AF mechanisms, with a combination of anatomical reentries, rotors, and meandering wavelets, but do not consider every proposed AF mechanism. For instance, focal activity was not simulated in our data, as it tends to terminate if its source is ablated. For example, the pulmonary veins, a common source of focal activity, especially in paroxysmal AF, would likely still require dedicated ablation even if our virtual ablation lines were clinically feasible.
The epi-endo dissociation mechanism was also not examined, since our work focuses on modeling the atrium as a surface of negligible depth. That mechanism requires each atrium to be treated as at least two surfaces with communication channels between them, or possibly even as a more complex 3D structure with volume.
Finally, we do not consider bi-atrial arrhythmias, which have similarities with the epi-endo case, as the two atria can be modeled as two surfaces with connecting channels between them. We are currently developing a new study that focuses on that topic.
5 Conclusion
This work demonstrated that much like in AT, the Index Theorem must also hold for simulated AF, using a large and realistic dataset of patient-based simulations, which brings more insight to the mechanisms of AF. We also showed that virtual ablation that simply connects opposite-chirality singularities may not terminate AF, highlighting issues with current ablation strategies. The greater aim was to guide future clinical work, dispelling misconceptions about functional reentries and guiding future ablation methods.
Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.
Author contributions
AB: Conceptualization, Formal Analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review and editing. RV: Conceptualization, Data curation, Investigation, Methodology, Software, Validation, Writing – review and editing. BV: Conceptualization, Methodology, Software, Visualization, Writing – review and editing. SL: Methodology, Software, Visualization, Writing – review and editing. AO: Methodology, Software, Visualization, Writing – review and editing. TN: Methodology, Software, Visualization, Writing – review and editing. SH: Methodology, Software, Visualization, Writing – review and editing. VS: Funding acquisition, Project administration, Resources, Writing – review and editing. NV: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review and editing.
Funding
The authors declare that financial support was received for the research and/or publication of this article. This research was supported by VLIR/iBOF project “DIAMOND” (Grant 20-VLIR-iBOF-027), and a “Starting Grant” from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant No. 900008), awarded to Nele Vandersickel. Vincent Segers was funded by the Senior Clinical Investigator fellowship (Application numbers 1842219N, 1842224N).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The handling editor OA declared a past co-authorship with the author NV.
Generative AI statement
The authors declare that Generative AI was used in the creation of this manuscript. During the preparation of this work, the authors used ChatGPT (GPT-4o) in order to enhance the readability of this work. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
Abeele R. V. D., Hendrickx S., Carlier N., Wülfers E. M., Santos Bezerra A., Verstraeten B., et al. (2025). DGM-TOP: automatic identification of the critical boundaries in atrial tachycardia. Front. Physiology 16, 1563807. doi:10.3389/fphys.2025.1563807
Allessie M., de Groot N. (2014). CrossTalk opposing view: rotors have not been demonstrated to be the drivers of atrial fibrillation. J. physiology 592 (15), 3167–3170. doi:10.1113/jphysiol.2014.271809
Arno L. (2021). “Phase defect lines during cardiac arrhythmias: from theory to experiment,” in arXiv preprint arXiv:2101.00315.
Arno L., Kabus D., Dierckx H. (2024). Analysis of complex excitation patterns using Feynman-like diagrams. Sci. Rep. 14, 28962. doi:10.1038/s41598-024-73544-z
Atienza F., Climent A. M., Guillem M. S., Berenfeld O. (2015). Frontiers in non-invasive cardiac mapping: rotors in atrial fibrillation-body surface frequency-phase mapping. Card. Electrophysiol. Clin. 7 (1), 59–69. doi:10.1016/j.ccep.2014.11.002
Courtemanche M., Ramirez R. J., Stanley N. (1998). Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am. J. Physiology-Heart Circulatory Physiology 275 (1), H301–H321. doi:10.1152/ajpheart.1998.275.1.h301
DeTal N., Kaboudian A., Fenton F. H. (2022). Terminating spiral waves with a single designed stimulus: teleportation as the mechanism for defibrillation. Proc. Natl. Acad. Sci. 119, e2117568119. doi:10.1073/pnas.2117568119
Duytschaever M., Van den Abeele R., Carlier N., Bezerra A. S., Verstraeten B., Lootens S., et al. (2024). Atrial topology for a unified understanding of typical and atypical flutter. Circulation Arrhythmia Electrophysiol. 17, e013102. doi:10.1161/CIRCEP.124.013102
Duytschaever M., De Smet M., Martens J., El Haddad M., De Becker B., Francois C., et al. (2025). How a topological mindset may offer extra control during mapping and ablation of left-sided reentrant atrial tachycardia. Circulation Arrhythmia Electrophysiol. 18 (7), e013780. doi:10.1161/CIRCEP.125.013780
Gurevich D. R., Grigoriev R. O. (2019). Robust approach for rotor mapping in cardiac tissue. Chaos An Interdiscip. J. Nonlinear Sci. 29, 5. doi:10.1063/1.5086936
Haïssaguerre M., Sanders P., Hocini M., Hsu L. F., Shah D. C., Scavée C., et al. (2004). Changes in atrial fibrillation cycle length and inducibility during catheter ablation and their relation to outcome. Circulation 109 (24), 3007–3013. doi:10.1161/01.CIR.0000130645.95357.97
Haissaguerre M., Hocini M., Denis A., Shah A. J., Komatsu Y., Yamashita S., et al. (2014). Driver domains in persistent atrial fibrillation. Circulation 130 (7), 530–538. doi:10.1161/CIRCULATIONAHA.113.005421
Lee S., Sahadevan J., Khrestian C. M., Cakulev I., Markowitz A., Waldo A. L. (2015). Simultaneous biatrial high-density (510–512 electrodes) epicardial mapping of persistent and long-standing persistent atrial fibrillation in patients: new insights into the mechanism of its maintenance. Circulation 132 (22), 2108–2117. doi:10.1161/CIRCULATIONAHA.115.017007
Lin Y.-J., Lo M. T., Lin C., Chang S. L., Lo L. W., Hu Y. F., et al. (2013). Prevalence, characteristics, mapping, and catheter ablation of potential rotors in nonparoxysmal atrial fibrillation. Circulation Arrhythmia Electrophysiol. 6, 851–858. doi:10.1161/CIRCEP.113.000318
Lootens S., Janssens I., Van Den Abeele R., Wülfers E. M., Bezerra A. S., Verstraeten B., et al. (2024). Directed graph mapping exceeds phase mapping for the detection of simulated 2D meandering rotors in fibrotic tissue with added noise. Comput. Biol. Med. 171, 108138. doi:10.1016/j.compbiomed.2024.108138
Marcotte C. D., Grigoriev R. O. (2017). Dynamical mechanism of atrial fibrillation: a topological approach. Chaos An Interdiscip. J. Nonlinear Sci. 27, 093936. doi:10.1063/1.5003259
Maury P., Takigawa M., Capellino S., Rollin A., Roux J. R., Mondoly P., et al. (2019). Atrial tachycardia with atrial Activation Duration Exceedang the Tadhycardia Cycle Length: Mechenisms and Prevalente. cACC Clin. Electrophysiol. 5 (8), 907–916. doi:10.1016/j.jacep.2019.04.015
Morgan R., Colman M. A., Chubb H., Seemann G., Aslanidi O. V. (2016). Slow conduction in the border zones of patchy fibrosis stabilizes the drivers for atrial fibrillation: insights from multi-scale human atrial modeling. Front. Physiology 7, 474. doi:10.3389/fphys.2016.00474
Nagy S. Z., Kasi P., Afonso V. X., Bird N., Pederson B., Mann I. E., et al. (2021). Cycle Length Evaluation in Persistent Atrial Fibrillation Using Kernel Density Estimation to Identify Transient and Stable Rapid Atrial Activity. Cardiovasc. Eng. Technol. 13 (2), 219–233. doi:10.1007/s13239-021-00568-1
Narayan S. M., Baykaner T., Clopton P., Schricker A., Lalani G. G., Krummen D. E., et al. (2014). Ablation of rotor and focal sources reduces late recurrence of atrial fibrillation compared with trigger ablation alone: extended follow-up of the CONFIRM trial (Conventional Ablation for Atrial Fibrillation With or Without Focal Impulse and Rotor Modulation). J. Am. Coll. Cardiol. 63 (17), 1761–1768. doi:10.1016/j.jacc.2014.02.543
Nattel S., Xiong F., Aguilar M. (2017). Demystifying rotors and their place in clinical translation of atrial fibrillation mechanisms. Nat. Rev. Cardiol. 14, 509–520. doi:10.1038/nrcardio.2017.37
Plank G., Loewe A., Neic A., Augustin C., Huang Y. L., Gsell M. A. F., et al. (2021). The openCARP Simulation Environment for Cardiac Electrophysiology. Under Rev. bioRxiv Prepr. 208, 106223. doi:10.1016/j.cmpb.2021.106223
Rappel W.-J., Baykaner T., Zaman J., Ganesan P., Rogers A. J., Narayan S. M. (2024). Spatially Conserved Spiral Wave Activity During Human Atrial Fibrillation. Circulation Arrhythmia Electrophysiol. 17, e012041. doi:10.1161/CIRCEP.123.012041
Roney C. H., Sim I., Yu J., Beach M., Mehta A., Alonso Solis-Lemus J., et al. (2022). Predicting Atrial Fibrillation Recurrence by Combining Population Data and Virtual Cohorts of Patient-Specific Left Atrial Models. Circulation Arrhythmia Electrophysiol. 15.2, e010253. doi:10.1161/CIRCEP.121.010253
Santucci P. A., Bhirud A., Vasaiwala S. C., Wilber D. J., Green A. (2024). Identification of 2 Distinct Boundaries Distinguishes Critical From Noncritical Isthmuses in Ablating Atypical Atrial Flutter. Clin. Electrophysiol. 10.2, 251–261. doi:10.1016/j.jacep.2023.09.024
Tomii N., Yamazaki M., Ashihara T., Nakazawa K., Shibata N., Honjo H., et al. (2021). Spatial phase discontinuity at the center of moving cardiac spiral waves. Comput. Biol. Med. 130, 104217. doi:10.1016/j.compbiomed.2021.104217
Vandersickel N. (2024). Impact of topology on the number of loops during macro-re-entrant atrial tachycardia.
Vandersickel N., Hendrickx S., Van Den Abeele R., Wuelfers E., Bezerra A. S., Fuenmayor S., et al. (2023). Unique topological classification of complex reentrant atrial tachycardias enables optimal ablation strategy. Europace 25, euad122.078. doi:10.1093/europace/euad122.078
Vigmond E., Pashaei A., Amraoui S., Cochet H., Hassaguerre M. (2016). Percolation as a mechanism to explain atrial fractionated electrograms and reentry in a fibrosis model based on imaging data. Heart Rhythm 13 (7), 1536–1543. doi:10.1016/j.hrthm.2016.03.019
Waks J. W., Josephson M. E. (2014). Mechanisms of Atrial Fibrillation – Reentry, Rotors and Reality. Arrhythmia and Electrophysiol. Rev. 3.2, 90–100. doi:10.15420/aer.2014.3.2.90
Wartenberg D. (1985). Multivariate Spatial Correlation: A Method for Exploratory Geographical Analysis. Geogr. Anal. 17, 263–283. doi:10.1111/j.1538-4632.1985.tb00849.x
Keywords: atrial fibrillation, cardiac electrophysiogy, reentrant arrhythmia, topology, index theorem
Citation: Bezerra AS, Van Den Abeele R, Verstraeten B, Lootens S, Okenov A, Nezlobinsky T, Hendrickx S, Segers VFM and Vandersickel N (2025) Reentry-driven model of atrial fibrillation is maintained by paired reentries and terminated by strategic pairwise virtual ablation. Front. Physiol. 16:1695431. doi: 10.3389/fphys.2025.1695431
Received: 29 August 2025; Accepted: 30 October 2025;
Published: 17 November 2025.
Edited by:
Oleg Aslanidi, King’s College London, United KingdomReviewed by:
Jason D. Bayer, Université de Bordeaux, FranceJosselin Duchateau, Centre Hospitalier Universitaire de Bordeaux, France
Copyright © 2025 Bezerra, Van Den Abeele, Verstraeten, Lootens, Okenov, Nezlobinsky, Hendrickx, Segers and Vandersickel. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Arthur S. Bezerra, YXJ0aHVyLmJlemVycmFAdWdlbnQuYmU=
Bjorn Verstraeten1