A Phase Defect Framework for the Analysis of Cardiac Arrhythmia Patterns

During cardiac arrhythmias, dynamical patterns of electrical activation form and evolve, which are of interest to understand and cure heart rhythm disorders. The analysis of these patterns is commonly performed by calculating the local activation phase and searching for phase singularities (PSs), i.e., points around which all phases are present. Here we propose an alternative framework, which focuses on phase defect lines (PDLs) and surfaces (PDSs) as more general mechanisms, which include PSs as a specific case. The proposed framework enables two conceptual unifications: between the local activation time and phase description, and between conduction block lines and the central regions of linear-core rotors. A simple PDL detection method is proposed and applied to data from simulations and optical mapping experiments. Our analysis of ventricular tachycardia in rabbit hearts (n = 6) shows that nearly all detected PSs were found on PDLs, but the PDLs had a significantly longer lifespan than the detected PSs. Since the proposed framework revisits basic building blocks of cardiac activation patterns, it can become a useful tool for further theory development and experimental analysis.


INTRODUCTION
About once per second, a wave of electrical depolarization travels through the heart, coordinating its mechanical contraction. The heart is a prime example of a dynamical system that can self-organize across different scales: changes in the flow of ions through the cell membrane affect the dynamics of the emergent pattern and may lead to life-threatening heart rhythm disorders. For more than a century, great efforts have been allocated to understand the different dynamical mechanisms of arrhythmia initiation and maintenance in the heart. It was found that electrical activity during cardiac arrhythmias may travel in a closed circuit within the cardiac wall, and thus re-excite the heart (Mines, 1913). Later on, Allessie et al. demonstrated with electrode recordings that re-excitation of the tissue by the wave itself was also possible without a central obstacle (Allessie et al., 1973). These rotating vortices, also called rotors, were later demonstrated using optical mapping techniques, both in ventricular tachycardia (Gray et al., 1995a) and ventricular fibrillation (Gray et al., 1998). Gray et al. (1998) also developed the concept of phase analysis to identify structures in the excitation patterns in the heart. They considered two fields V( r, t) and R( r, t) of the system and calculated the phase of a point of the medium as the polar angle in the (V, R)-plane, see Figure 1: (1) An additive constant in Equation (1) is included here to fix the absolute phase, e.g., such that resting state corresponds to φ act = 0. To discriminate from an alternative phase definition below, we call this classical phase the "activation phase" and denote it as φ act . The application of phase analysis to cardiac electrical signals during arrhythmia showed that there are few points in the medium where all possible phases meet each other, such that the point itself has no well-defined phase, see Gray et al., 1998 and Figure 1C. Such point is now generally accepted to be a phase singularity (PS), and since then PSs have been widely used in the analysis of cardiac excitation patterns. In two dimensions (2D), PSs are associated with the rotor core, i.e., the regions around which electrical waves revolve tachycardia or fibrillation.
In three dimensions (3D), the set of PSs in the medium extend to a dynamical curve, the rotor filament. The rotor filament has been used in many modeling and theoretical studies, as it allows for easy visualization of the rotor dynamics (Wellner et al., 2002;Clayton et al., 2005;Verschelde et al., 2007;Dierckx et al., 2012).
To identify PSs in datasets, different methods have been developed (Fenton and Karma, 1998;Gray et al., 1998;Bray and Wikswo, 2003;Clayton et al., 2005;Kuklik et al., 2015Kuklik et al., , 2017. However, these methods that look for PSs have limitations when performing detections near conduction block lines, as illustrated in Figure 2. A first limitation of the classical PS concept is non-robustness of PS detection methods. Figure 2A shows an example of a simulation of a 2D sheet of cardiac tissue using the S1-S2 stimulation protocol. The second (S2) pulse undergoes unidirectional block in the wave of the S1 pulse, leading to the formation of a conduction block line (CBL), rendered in black. The region excited by the S2 stimulus grows, turning around the CBL endpoints, but no PSs are detected at its end. Later on, when the initial S2 region has repolarized, a large number of PSs is found at the CBL, even though there is no rotor core. Moreover, FIGURE 1 | Classical phase analysis of cardiac rotors. (A) Circular-core rotor with Aliev-Panfilov kinetics (Aliev and Panfilov, 1996), where in each point an activation variable u and a recovery variable v are defined. (B) Two observables of the system, V = u and R = v plotted against each other reveal a cycle corresponding to the action potential. The polar angle with respect to a point (V * , R * ) situated within the cycle serves as a definition of activation phase. (C) Coloring the rotor using phase and a periodic colormap reveals a special point where all phases meet, the phase singularity (PS), see white dot. This point has V = V * , R = R * . the precise number of PSs found depends on the algorithm used, as shown in the third and fourth panel of Figure 2A. In literature, other false positives of PS identification have been reported (Podziemski et al., 2018).
A second limitation of classical PS detection algorithms is related to rotors with so-called linear cores, i.e., rotors which have a CBL in the central region (core) around which they rotate, see Figure 2B. If one looks at the point where the wave front ends on the CBL, this is in the classical picture a PS, see Figure 2B. However, looking closer, only 3 phases are found in the neighborhood of this point: excitable tissue (ahead of the front, blue), recently excited tissue (behind the front, green), and refractory tissue (across the CBL, red). A point touching three phases is not necessarily a mathematical PS, as the latter should touch all phases ( Figure 5).
If we represent the same linear core rotor in terms of the local activation time (LAT), as indicated in the rightmost panel of Figure 2B, an extended CBL is clearly seen, where isochrones of activation coincide. However, this extended CBL is not captured by the classical PS analysis.
A third limitation of the classical PS detection methods concerns the coincidence of PSs and CBLs. Podziemski et al. (2018) and our own optical mapping experiments of ventricular tachycardia in rabbits indicate that the PSs are generally localized at a CBL, see Figure 2C. When comparing different time frames, we moreover find that the PSs tend to jump between different locations on the CBL, while the CBL itself seems to persist longer in time.
The goal of this paper is to present a new framework for phase analysis, which will lead to identification of novel structures in the regions of conduction block: phase defect lines (PDLs) and phase defect surfaces (PDSs) instead of PSs and filaments. In the context of cardiac arrhythmia patterns, interpreting linear rotor cores as a phase defect was only done very recently by Tomii et al. (2021) and in our preprint to this manuscript (Arno et al., 2021). Within this paper, we additionally introduce a phase based on local activation time, which bridges the gap between the phase viewpoint and LAT viewpoint. We propose a simple PDL detection method and apply it to data from simulations FIGURE 2 | Limitations of current PS detection algorithms. (A) Application of the S1-S2 stimulation protocol in the BOCF (Bueno-Orovio et al., 2008) model to initiate a rotor. White dots denoted detected PSs with either a 2 × 2 ring (third panel) or 2 × 2 + 4 × 4 ring (fourth panel) using the method of Kuklik et al. (2017). Both methods identify multiple PS on the CBL (black line). (B) A linear-core rotor in the Fenton Karma model, from 3 perspectives: transmembrane voltage (left), activation phase (middle panel) with PS indicated (white) and LAT (right). (C) Optical mapping of rabbit hearts during ventricular tachycardia showing that detected PSs are all located on CBLs.
(in 2D and 3D) and optical mapping experiments. We draw conclusions from the analysis of our experiments and outline how future theory development can be initiated from the phase defect concept.

Numerical Simulations
To illustrate our theoretical concepts, numerical simulations were performed in a cardiac monodomain model: (2) This partial differential equation states how a column matrix u( r, t) of m state variables evolves in time. The constant matrix P = diag(P 11 , 0, ...0) makes sure that only the transmembrane potential u 1 undergoes diffusion. Note that our analysis methods are applicable to excitation patterns in general and are not restricted to reaction-diffusion systems only.
In this work, we use the reaction-kinetics functions F(u) from Aliev and Panfilov (1996) Fenton and Karma (1998) (FK model,m = 3), and Bueno-Orovio et al. (2008) (BOCF model,m = 4). The equations are numerically integrated using the Euler method with the spatial resolution dx and time step dt as specified in Table 1. Only in the example with a biventricular geometry, uniaxial anisotropy was included: the Laplacian in (2) was replaced by f e j f , with D 1 = 1, D 2 = 1/5 and e f the local fiber orientation as measured by Hren and Stroink (1995).
The first and second state variables in these different models were used to calculate φ act , see Table 1 for definitions and threshold values.  Aliev and Panfilov, 1996Fenton and Karma, 1998Bueno-Orovio et al., 2008 Parameter set Standard MLRI(2D) and MBR(3D) Epicardial (EPI)  ,5,6,8,11 Figures 2,4,7,9,12 FIGURE 3 | Method of PS detection by Kuklik et al. (2017) as used in this paper. Phase differences are considered between pairs of points on a ring of 4 points (Left, 2 × 2 method) or 8 points (Right, 4 × 4 method), to assess whether a PS is present at the central location.

Optical Mapping Experiments
We applied our new framework to experimentally recorded movies of rotors in isolated Langendorff-perfused rabbit hearts (n = 6), as described by Kulkarni and Tolkacheva (2018).
Optical movies corresponding to the fluorescence signal were recorded from the epicardial surfaces of the left or right ventricular surface of the heart (LV and RV) at 1,000 frames per second, with 14-bit, 80 x 80-pixel resolution cameras (Little Joe, RedShirt Imaging, SciMeasure) after a period of stabilization of approximately 30 min. Data processing was performed using a custom-made program in Matlab 1 . The background was removed by thresholding, and in each pixel the intensity was normalized against the minimal and maximal value of optical intensity attained in that pixel over the full recording. The characteristic timescale of activation was found as the inverse of 1 Matlab. version 9.9.0. (2020). The MathWorks Inc., Natick, Massachusetts.
the dominant frequency of activation was used to evaluate the Hilbert transform to obtain the activation phase φ act (Bray and Wikswo, 2002).

Methods for PS and CBL Detection
PSs were detected using the method of Kuklik et al. (2017), see Figure 3. A ring of N = 4 or N = 8 pixels was considered around each grid cell of the medium, and the phase difference was computed between adjacent points on the ring. If exactly one of these phase differences was larger than π in absolute value, a PS was assigned to that grid element: In the 2 × 2 method, only the 4 points around each grid cell are used; in the 2 × 2 + 4 × 4 method, a second ring of 8 pixels was used and a PS is assigned if both the small ring and large ring have one phase difference larger than π, see Figure 3. To localize CBLs, frames were recorded every t time units, and if in the new frame, V rose above V * , this time was locally saved as the newest LAT. If an edge connecting two neighboring points of the grid had two LAT that differed by more than t c ≈ dx/c (with c the plane wave velocity in the medium), the middle of that edge was considered to be part of a CBL: |t activation ( r 1 ) − t activation ( r 2 ))| > t c and t activation ( r 1 ) > 0 and t activation ( r 2 ) > 0 Here, the first condition selects the union of the WF and the CBL.
The threshold values t c that we used with the different model kinetics are given in Table 1. The second and third condition are imposed to obtain CBLs only.

Introduction of Phase Defects
In many physical sciences, interfaces are found between regions of different phase, where they are known as domain walls or phase defects. In this sub-section, we demonstrate the existence of phase defects. Figure 4 shows in each row the time evolution of a simulated linear-core rotor (at four different times). When applying classical phase analysis with φ act , see Figure 4B, we note that there are several interfaces that connect excited ( Figure 4A, yellow) to non-excited ( Figure 4A, blue) areas. These interfaces can be shown by plotting the spatial gradient of φ act , see Figure 4C. When taking the spatial gradient, it is important to disregard phase differences of 2π. We compute spatial derivatives of phase in the grid point with discrete position (i, j) as: with U the unwrap function, which adds an integer multiple of 2π to its argument, to bring the result closest to zero: In the linear-core rotor, we note different steep transitions between zones of approximately equal phase, and discriminate between them as follows. We denote by φ 1 a phase in the middle FIGURE 5 | Closer look at a linear-core rotor. At the location where classical methods detect a PS, three distinct phases come together: recovered, excited and refractory tissue. At either side of the CBL, two distinct phases are present: refractory vs. either excitable or recovered. Therefore, the CBL is a phase defect line (PDL).
of the upstroke (depolarization), and φ 2 a phase in the middle of the down stroke (repolarization). These values are found as Thereafter, the wave front (WF) and wave back (WB) can be defined as (with ♥ depicting the domain or cardiac tissue): In the classical view, the only way to connect WF and WB is in a PS, where different phases φ 1 and φ 2 can meet. However, if a wave hits unrecovered (refractory) tissue, unidirectional conduction block occurs and a CBL is formed. In that case, a third kind of interface is seen, in addition to WF and WB, see Figure 5. Since the tissue at either side of a CBL was activated at a different time, it also has different phase, see different colors across the black line in Figure 5. This observation implies that not only the end point of a wave front (the classical PS) is a special point, but that all points on a CBL have a nontrivial phase structure: they are situated on a region where phase suddenly transitions from one value to another, which we call a phase defect line (PDL).
Note that the PDL shown here persists in time until either both sides have fully recovered, or when one side of the line is re-excited again, while the other is not. Figure 4C shows the norm of ∇φ act , which visualizes the domain walls between the different zones of similar phase; note that without further filtering, the transition zones consist of the union of WF, WB, and PDL.

An Alternative Phase Definition Based on LAT
We find it useful to discriminate WF and WB from CBLs, since the former are "natural" behavior in an excitable system, while conduction block can lead to the formation of abnormal patterns and arrhythmias and therefore, deserves to be detected separately. Therefore, we now introduce another definition of phase, which is only sensitive to CBLs.
As we described before, the classical definition of phase is not unique, as one could make various choices for the observables V and R, and their threshold values V * , R * can be also chosen at will, as long as they are located within the cycle in the (V, R)-plane.
The most important feature of the phase is its periodicity. For convenience, we choose the phase of the resting state to be phase equal to 0. Apart from that, there is need to let the phase correspond to the polar angle in one choice of (V, R)-coordinates: we can define new equivalent phases φ by a continuous mapping of φ act : The last condition ensures that the transformation between the phases is invertible. In a mathematical sense, fixing a phase is nothing but introducing coordinates on the inertial manifold (dynamical attractor) on which the dynamics takes place. As with coordinates in the plane or on a surface, many choices are possible, and depending on the circumstances, some may be more appropriate than others.
To link the concept of LAT with phase, we propose to take the elapsed time as the new phase, since in leading order (neglecting inhomogeneity, electrotonic effects and long-term memory), the state of a point will depend on the elapsed time since previous activation. However, t elapsed ∈ R + , which is an unbound interval, such that t elapsed itself is not a suitable replacement for φ act . To this purpose, we apply a sigmoidal function to t elapsed , in order to define an "arrival time phase": φ arr ( r, t) = 2π tanh(3t elapsed ( r, t)/τ ).
Here, τ is a characteristic time constant of the medium, e.g., the mean action potential duration. Figures 6A,B show how φ arr depends on both t and φ act for FK reaction kinetics. From this, it can be seen that, to a good approximation, φ arr is a reparameterization of φ act . The new phase φ arr is not constant in time, but gets update at the places where the wave front excites tissue, just like the updates in LAT. It relies on two parameters, which we chose manually for now: the timescale τ and the critical threshold V * to consider tissue as excited, see definition (8). Importantly, φ arr unifies the classical phase description (since φ act is a phase) with the LAT description (since it is defined based on LAT). When a linear-core rotor is plotted using φ arr , as shown in Figure 7, we see that the WF and WB are not strongly emphasized anymore, since both occur near φ arr = 0 (omitting phase differences of 2π). Still visible as a spatial phase transition is the region of conduction block, or the PDL, see Figure 6C.
The "wave arrival time phase" φ arr is thus a kind of hybrid between classical phase and the LAT viewpoint. For, the lines of equal phase will now correspond to the isochrones that commonly appear in medical literature and our Figure 2B. At  Figure 2B, now shown with φ arr . Note that WF and WB are no longer showing abrupt phase variations, these only happen at the PDL, i.e., the points where conduction block happened. (D) Magnitude of the gradient of φ arr . Due to the discrete sampling of LAT, a staircase artefact in the gradient is seen (no smoothing was applied here). a CBL, LAT changes discontinuously across the conduction block site, and this phase defect region now clearly exhibits a discontinuous phase. So, depending on the definition of phase used (φ act or φ arr or another one) the phase defect will be either a steep but continuous phase transition, or a true discontinuity in phase.

Interpretation of Phase Defects as Branch Cuts From Complex Analysis
Here, we demonstrate how phase defects can be interpreted as a branch cut from complex analysis (Arfken and Weber, 1995). For any non-zero complex number of the form z = x + iy with i 2 = −1, the phase can be calculated using Next, one can consider functions of a complex variable: w = f (z) with both w and z complex numbers, e.g., w = z 2 . To make a visual representation of the function value, one can draw colormaps of the phase arg(w), or show arg(w) as an XY-dependent elevation in a so-called Riemann surface; see Figure 8. In case of polynomial or rational functions, such phase map will reveal point singularities similarly to PSs found in circular-core spiral waves, see Figure 8A. When the phase of the simple function w = f (z) = z is shown as a graph of (x, y) the phase surface resembles a helicoid. If one walks around the PS once in the XY-plane, the phase will change by 2π, but since phase is measured disregarding terms of 2π, one ends in the same state. This effect is more easily seen using a cyclic colormap in 2D, see Figure 8A. Close to the central axis of the helicoid, the height (phase) is undefined, which is typical for a PS. In contrast to polynomial and rational functions, functions f which include roots, exhibit a discontinuity in their phase, also called "branch cut" in mathematics. In particular, for the function w = f (z) = √ z 2 − 1, there is a line of discontinuous phase for (x, y) between (−1, 0) and (1, 0). This line is easily noticed in the representation as a color map as a sudden transition in color, see Figure 8. This situation is in our opinion very similar to the phase discontinuity in the linear-core rotor throughout this paper.
For, if one regards two neighboring points on both sides of a conduction block line at the rotor core, they will have distinct phases. There need not be a single point where all phases spatially converge (that is, a PS), as on both sides the phases can gradually change to take the same value at the points where the branch cut stops.
From this paragraph and Figure 8 it is apparent that the PSs and PDLs are different topological structures. We argue in this paper that we and Tomii et al. (2021) are the first to note this difference in a cardiac electrophysiology context, and that discriminating between these different structures can be useful in theory development and applications.  Rigidly rotating spirals, as in the Aliev-Panfilov (AP) reaction-diffusion model (Aliev and Panfilov, 1996) correspond to a PS (gray), similar to the mathematical function φ(x, y) = arg(z) shown to the right of it. (B) Linear-core cardiac models, e.g., (Fenton and Karma, 1998) exhibit a PDL or branch cut (black/gray), like the mathematical function w = f(z) = arg( √ z 2 − 1) shown to the right of it. Gray areas denote a jump in the phase over a quantity not equal to an integer multiple of 2π, i.e., a PDL (physics) or branch cut (mathematics).

Methods for Phase Defect Detection
The numerical identification of PDLs in cardiac excitation patterns can be performed using various methods. Based on a phase (either φ act or φ arr ), one can assign a PDL to the midpoint of edges along which the phase gradient exceeds a predefined value, taking into account that phase differences of 2π are excluded: However, if the LAT is known, this method approximates the more direct method of identifying a CBL, see Equation (4). A FIGURE 9 | Interpretation of a CBL as a phase defect clarifies why a PS detection finds PSs on it. (A) φ act , rendered in 2D and 3D as a Riemannian surface. (B) φ arr , rendered in 2D and 3D as a Riemannian surface. At the CBL/PDL, the phase surface has a cliff-like appearance.
FIGURE 10 | Creation of a PS from a PDL using a S1-S2 protocol with AP kinetics, showing coexistence of PDL and PS in models that generate circular-core spirals. Snapshots (A-D) shown at times t = 20 ms, t = 38 ms, t = 57 ms, and t = 150 ms.
comparison between different manners to numerically compute PDLs will be deferred to another publication. Below, we also compare the lifetime of PSs and PDLs in experimental recordings. To estimate PS lifetime, we calculated PSs at each time frame between 501 and 1,500 ms using the Hilbert transform of V(t), followed by PS detection using 2 × 2 + 4 × 4-method, which is considered as robust in literature (Kuklik et al., 2017). If a PS moved between two subsequent frames at most 1 pixel in the horizontal, vertical, or diagonal direction, while keeping its chirality, it was considered PS. Otherwise the PS was considered a new one.
A similar approach was taken for PDLs. PDLs in experiments were computed using Equation (13), with manually chosen threshold φ crit = 2.22. If the middle of an edge was on a PDL, both neighboring nodes of the Cartesian grid were considered to be part of the PDL. A PDL was assumed to persist in time if the PDL at the new time and at the old time had points which were not further than 1 pixel away from each other (in horizontal, vertical, or diagonal direction). PDLs were allowed to branch and merge, and the lifetime was computed as the earliest and latest time which were connected by this family of PDLs.
A PS was furthermore considered to lie on a PDL if its distance was less than 0.5 pixels from a point on the PDL. PDLs in Figures 2C, 13 were visually rendered using splines. The wave front was added to these figures using Equations (7)-(8), with V * = 0.5.

RESULTS
Here we illustrate the appearance and relation between PSs and PDLs in simulations and an optical mapping experiment.

Rendering of Cardiac Activation Phase Around a Conduction Block Line as Riemannian Surfaces
The CBLs presented in Figure 2 are examples of PDLs: to this purpose, we redraw the rightmost panel of Figure 2A here using φ act and φ arr , rendered in 2D and 3D in Figure 9. Here, it can be seen that the suspected PSs in fact lie on the line where phase is nearly discontinuous.
We distinguish two cases here. If no rotor is attached to the phase defect, the Riemann surface looks like a sheet of paper with a cut in it and both edges shifted relative to each other, see Figure 9. If a rotor is attached to the phase defect, the Riemann surface looks like two rising slopes connecting two floors in a parking lot, see lower right panel in Figure 8B.

Coexistence of PS and PDL in the Same System
Here we demonstrate that PDLs can also exist in media that sustain spiral waves with PSs. This is illustrated in simulations where application of S1S2 pacing induces a PDL at the WB of the first stimulus, see Figure 10. Since AP reaction kinetics were used (Aliev and Panfilov, 1996), the PDL will disappear after both sides have reached full recovery, resulting in a PS at the core of the spiral wave, see Figure 10D. Hence, we conclude that PDLs and PSs can coexist in media supporting circular-core spiral waves, since the PDL is found at conduction block sites and PSs are found in the spiral wave core.

Absence of a PS at the Edge of a PDL
Although Figure 10 shows a PS situated near the end of a PDL, this is certainly not always the case. Counter examples FIGURE 11 | Analysis of a break-up pattern in a 3D slab geometry with linear-core rotors and rotational anistropy in the FK model. Panels (A-C) show activity at times t = 30, 35, and 40 ms, respectively. The rows show respectively 3D activity with filaments (colors indicate different filaments), normalized transmembrane voltage on the bottom surface, and analysis of the surface patterns in terms of φ arr and φ act and 3D activity with filaments colored as above and PDSs indicated in gray. Note that several phase defects are observed, indicating conduction block near the spiral wave core.
of PSs located away from end points of a PDL can be seen in Figures 2, 4.
The underlying mathematical reason is the following. At the edge of a PDL, the phase surface starts to be "torn, " see Figure 8B. Still, at this point itself, the phase is well-defined. This situation is different from a PS, see Figure 8A, where the Riemann surface locally resembles a staircase surface and the phase itself cannot be determined.

Phase Defect Analysis of Scroll Wave Turbulence
The phase defect framework was illustrated above on simple examples, but is designed to analyse much more (complex) dynamics in excitation patterns. We conducted a 3D simulation in the FK model with rotational anisotropy-induced break-up, as detailed in Fenton and Karma (1998). Figure 11 shows snapshots of the turbulent pattern on the bottom surface of the slab (representing endo-or epicardial view in this geometrical model).
Here, it can be seen that multiple PSs are detected on CBLs, that correspond to the core of cardiac rotors.

Numerical Results in Three Dimensions
In classical 3D theory, the PS extends to a dynamically varying curve called the rotor filament. With linear cores in the phase defect framework, the PDL will also extend in three dimensions, to form a phase defect surface (PDS). The bottom row of Figure 11 shows the PDSs for the full 3D break-up patterns in the FK model. Here, it is seen that the PDLs that we identified on the myocardial surface extend to surfaces in three dimensions. From the surface recording (2D) already, it can be seen that these surfaces can branch and merge in time.
A second example of PDSs is shown in Figure 12A, where a scroll wave with linear core was initiated by setting initial conditions with unidirectional block in the left-ventricular and right-ventricular free wall and the intraventricular septum. When viewed in terms of phase defects, 3 PDSs are seen. The two PDSs on the left correspond to filaments, and the rightmost PDS is a CBL; no matching filament is seen in Figure 12A.

Analysis of Rotors From Isolated Rabbit Hearts
Electrical activity of a Langendorff-perfused rabbit heart was visualized during ventricular tachycardia via optical mapping experiments (Kulkarni and Tolkacheva, 2018) as detailed in section 2.1.2. The processed movies of the different experiments performed are accessible in the Supplementary Material, and various characteristics of rotors in these movies are summarized in Table 2.
A first important observation in Figure 13 and Table 2 is that in more than 99% of the cases where a PS was found, it was situated on a PDL. A typical sequence of the activation is presented in Figure 13. Here, a classical rotor can be recognized from the progression in time of a wave front (gray line) that starts in the top-right corner and circles counterclockwise around a line of conduction block that comprises the center of this rotor. Note that the PS detection algorithm finds a variable number of PSs on a PDL over time and therefore, the lifetime of the PSs are smaller than the mean PDL lifetime, as is shown in Table 2. In contrast, the mean lifetime of PSs was about 3 frames.
Second, two types of PDLs were observed: approximately two thirds of them were not associated with any PS during their lifetime and we recognize these as conduction block lines without a rotor attached to them. About one third of PDLs/CBLs had at least one PS on them during their lifetime, and these are therefore, at the rotor core.
A third observation is that no rotors in the observed ventricular tachycardia completed a full turn on the epicardial surface, confirming other results in literature; a representative example is shown in the rightmost panel of Figure 13. Prior to completing a full rotation, excitation is surfacing from deeper layers, leading to a full conduction block. After that, the excitation loops in the same sequence, leading to a nearly Here, PSs were computed using the integral method using the "2 × 2 + 4 × 4" ring of points (Kuklik et al., 2017). Phase was computed for both PSs and PDLs with V(t) the normalized optical intensity and R(t) its Hilbert transform.  periodic ventricular tachycardia. Note that the fact that the rotor is not performing a full rotation on the epicardial surface is not reflected by the PS analysis, but becomes visible only in the PDL framework. We verified visually that this scenario was true for most re-entrant activity in the data. Figure 14 shows that in every recorded frame the number of PSs correlates with the total length of PDLs in that frame (data from one experiment is shown, counting PSs observed with two cameras on both sides of the heart). The intersection of the trend line with the vertical axis at finite distance shows that even without PSs, PDLs are found, which we attribute to CBLs without a rotor connected to it.

DISCUSSION
In this manuscript, we proposed a topological framework for excitable systems that feature conduction block lines. We brought terminology from complex analysis (Riemann surfaces, branch cuts) to a cardiac electrophysiology context and demonstrated that the phase defect concept can describe more structures than the classical PS. Moreover, we introduced the arrival time phase, which allows to convert activation time measurements to phase, and thereby unifies the isochrone and phase description of cardiac tissue.
Our main findings are the following (1) the distribution of phase in our cardiac simulations and experiments is organized in regions with slow spatial variations of phase, separated from each other by localized interfaces, which we called phase defects; (2) that when arrival time phase is used, the only aberrant interfaces occur at CBLs; (3) at these PDLs, classical detection methods tend to localize PSs.

PDLs vs. PSs
The concept of extended PDLs can in our opinion explain some limitations associated with classical PS analysis that we presented in Figure 2. First, near the central region of a PDL, the phase changes abruptly, and if this phase difference exceeds π, these points on a PDL can be classified by traditional methods as PSs, see Figure 2A. One option is to filter the resulting PSs and remove closely spaced PSs, but another option suggested here is to recognize the structures as line defects rather than point defects of phase.
The observation that points on a PDL can be detected as PSs at once explains why PSs are often found on CBLs, as was seen in Figure 2C. Here, we conclude that in systems that form linearcore rotors (even short-lived) do not possess PSs, but only PDLs. Conversely, systems that sustain circular-core rotors can have PSs at the spiral core, and PDLs at conduction block lines.
In Figure 2B, we argued that a PS cannot capture the shape of a linear core. In contrast, a PDL can do so, as it is a line spanning the entire conduction block region.
Geometrically speaking, PDLs and PSs are related to each other: a PS can be seen as the limit of a PDL having length zero. Conversely, a PDL can also be regarded as some kind of extended PS. Hence, the PDL concept can be seen as a refinement of a PS. Therefore, we make the suggestion to adapt current phase analysis methods to accommodate for PDLs, which can potentially make analysis more robust and accurate.

Physiological Interpretation of Phase Defects
In analogy to the term "phase singularity" referring to spiral wave tips, we here suggested the term "phase defect line" (PDL) to refer to CBLs when their phase structure needs to be emphasized. PDLs are ubiquitous in other branches of physics, e.g., magnetism (Landau and Lifshitz, 1935), liquid crystals (Williams, 1963), and string theory (Vilenkin, 1985).
Our justification as to why PDLs exist in these systems is based on biological rather than mathematical arguments and essentially traces back to the arguments of Winfree (1974). Consider the activation cycle of an excitable cell, i.e., its action potential, as shown in Figure 1B. Now, if one expects a PS in the center of a vortex, one assumes that the cells can also be in a state that lies somewhere in the middle of the cycle in state space. However, this situation may be biologically impossible: electrophysiology processes during activation will push the cell along its activation cycle, not necessarily allowing it to occupy the middle state. In case of a PDL, the "forbidden" state does not need to be realized. Going back from state space to real space, the region with forbidden states is mapped to the PDL.
If one, however, models cardiac tissue in the continuum approximation using a diffusion term, the transmembrane potential becomes continuous, and if this potential is used as one of the observables (V, R) to infer phase, points in the "forbidden zone" of state space are actually realized by the smoothing effect, leading to a finite thickness of border between regions of different phase. Since this diffusive effect acts during one action potential with duration τ , we estimate the effective thickness of the phase defect line to be with D the diffusion coefficient for transmembrane potential in the direction perpendicular to the PDL. For the BOCF model for human ventricle, one has D = 1.171 cm 2 /s and τ = 269 ms, leading to a PDL thickness of d = 5.6 mm if the PDL is perpendicular to the myofiber direction. With conduction velocity anisotropy ratio c 1 /c 2 ≈ 2.5, one finds d = 14 mm if the PDL is oriented parallel to the myofiber direction. These are rough estimates in a minimal monodomain model. Thus, as with other interfaces between phases in nature, a boundary layer forms with a thickness related to fundamental constants of the problem. This boundary layer effect explains why the common PS detection methods usually do return a result, even in the region of conduction block. When the new framework is applied to a linear-core rotor, the classical spiral wave tip or PS is expected to be the point where a wave front ends on a PDL. However, this may depend on the precise tip structure, spatial resolution and PS detection method, see e.g., Figure 4.

Relation to Other Works
Between the submission of our preprint (Arno et al., 2021) and acceptance of this paper, a paper was published by Tomii et al. (2021), who also identify a phase discontinuity at the center of cardiac vortices, and compares it to a branch cut in complex analysis. This effect was demonstrated by Tomii et al. in numerical simulations only. They detected PDLs using This method is very similar to our Equation (13). However, our study goes further than their result, as we introduce arrival time phase to link the phase picture with LAT, we demonstrate the presence of PDLs in optical mapping experiments, and report about PDSs in 3D as well.
We are also aware that the concept itself of a CBL is not new at all, as conduction blocks are mentioned for more than a century in electrophysiology reports (Mines, 1913). In 3D, filaments have also been reported as being ribbon-like before (Efimov et al., 1999). However, these observations have not been thoroughly integrated in pattern analysis methods (which is usually based on PSs) or the dynamical systems approach to better understand arrhythmias (which has long considered circular-cores only).
Meanwhile, CBLs are routinely localized and visualized during clinical procedures, e.g., as the result of an ablation. A prime insight advocated in this manuscript, and by Tomii et al. (2021), is that the CBL itself is an extended line of abrupt change between two regions with different activation phase, and that this resembles more a mathematical branch cut than a PS. Further studies are recommended to see if the phase defect concept can help to solve some basic questions in the field, such as why rotors do not complete full turns in tissue, and which characteristics of the reaction kinetics determine the core shape.

Relevance of Phase Defects for Theory of Arrhythmias
The recognition that a PS may not be the best option to describe linear core rotors, is in our opinion opening up several paths for further analysis and insight in cardiac arrhythmia patterns.
First, by introducing the arrival time phase φ arr , we have unified the LAT description with the phase description, which allows computing φ act . Which phase to use in practice will depend on the availability of data and their quality. Figure 6 furthermore shows that both phases can be converted into each other. Hence, activation phase can also be estimated from LAT, and arrival time phase from a single snapshot of V and R at a given time.
Second, this work was inspired by localizing the regions where an external stimulus can affect rotor dynamics. The relevant sensitivity function is here known as a "response function" (RF) (Biktasheva and Biktashev, 2003), and these have been computed recently also for meandering and linear-core rotors (Marcotte and Grigoriev, 2015;Dierckx et al., 2017). For circular-core rotors, it has been numerically demonstrated that the RF is located near the spiral wave tip, or its PS. So, rather than acting as a point particle, the rotor has a finite extent given by the spatial decay width of its RFs, with sensitivity localized near the spiral wave tip or PS, not at its instantaneous rotation center. For linear-core spirals, the PDL at their core indicates that the sensitivity is located around this zone of rapidly varying phase. More specifically, their RF was shown to be localized near the end points of the PDL (Dierckx et al., 2017) ("turning points"), suggesting that linear-core rotors resemble a localized particle that hops from turning point to turning point. That the sensitivity is lower in the middle part of the PDL can be understood from the conduction block interpretation: this block line cannot move much, as the tissue on one side is inexcitable. Hence, the PDL concept confirms insights from RF theory, and it can be hoped that combining both may open up ways to deepen the understanding of how point dynamics (local electrophysiology) affect the emerging patterns.
Third, it is yet unclear how essential structures such as WBs and CBLs move under the influence of external stimuli such as impeding waves or electrotonic currents from boundaries or 3D effects. To answer this question is a key step toward designing better methods to control wave patterns, e.g., as in defibrillation. The identification of the phase defect, with boundary layer, could allow in the near future to design a quantitative theory of their dynamics, in analogy to earlier successes on rotor filaments with circular cores (Keener, 1988) and excitation fronts (Kuramoto, 1980), which have lead to instability criteria for pattern formation in the heart (Biktashev et al., 1994;Dierckx et al., 2012).
Fourth, a consequence of this investigation is that describing fibrillation patterns using PSs only (Gray et al., 1995b;Gurevich and Grigoriev, 2019) may not be sufficient to fully understand the dynamics if conduction block occurs. Therefore, these studies could benefit from additional identification of PDLs, of which we made first steps in Figures 11, 13.
Finally, our preliminary results on 3D dynamics in Figure 11, and Figure 12 is a first step to extending the filament concept to the PDSs, like was done here to go from PSs to PDLs. In 2D, we here demonstrated that circular-core spiral waves are centered around a PS and linear-core meandering spirals rotate around a PDL. Since in 3D, the collection of classical PSs forms a filament curve, linear-core scroll waves extend to PDSs. Therefore, the study of the dynamics of meandering scroll waves can benefit from its description in terms of phase defect surfaces. These surfaces, as shown in Figures 11, 12 have a ribbon shape (Efimov et al., 1999), but the meaning of the ribbon is different from the ribbon model for filaments (Echebarria et al., 2006), where the surface normal of the ribbon was purely related to scroll wave twist. In our framework, the width of the PDS (ribbon) is approximately equal to the length of the PDL in 2D. How dynamical concepts such as filament tension (Biktashev et al., 1994) and rigidity (Dierckx et al., 2012) (Barkley et al., 1990), between the extremes of circular and linear cores, it still needs to be determined whether the phase structure is a PS, PDL or a different structure. Further analysis of complex 3D dynamics in the phase defect framework will be the topic of further investigation.

Interpretation of Experimental Results
Our analysis of rotors observed in isolated rabbit hearts during experiments shows in short two observations. First, that PSs on the epicardium were in this system always located on conduction block lines (PDLs). This observation is in line with older and more recent literature (Efimov et al., 1999;Podziemski et al., 2018). Second, the rotors were not performing more than one full rotation. To our knowledge, the question on why rotors are not performing multiple rotations in ventricular tissue is unresolved. The phase defect framework may be used to answer this question in a comprehensive manner, rather than resorting to numerical simulations only.
Another relevant question is whether the PS or PDL is the structure that governs the organization of cardiac arrhythmias. Here, we would answer that both are intricately related to each other: in systems with linear-core rotors, we think the PS only arises due to smoothing of the phase field and is therefore located at a specific spot on the PDL, often where it meets the wave front. Whether one should localize PSs or PDLs in experiments or clinical situation, depends on the application and desired accuracy.

Limitations and Outlook
The concept of PDLs was illustrated here in a very simple setting: we mostly considered single rotors or rotor creation in 2D. Our initial results in 3D (see Figure 12) and with breakup patterns (see Figure 11) show that there is rich dynamics present that can be further analyzed using the phase defect framework. We also did not consider alternans here, as in that case the inertial manifold in state space will be not a simple closed curve, and may require two or more phase parameters to describe those states. This effort will be undertaken in future work. Also, we have described the framework in a general setting, without explicit attention to the reason why rotors form in the system. Future work could be directed to elucidate differences between rotors formed by dynamical breakup (e.g., in atrial fibrillation models) or by interaction with inhomogeneities of the medium.
In this paper, the phase defect concept was used only to describe the observed patterns, without yet looking into how the topology determines the evolution of the structures such as WFs, WBs, and PDLs. Still, we developed the phase defect framework with the aim of providing a comprehensive quantitative analysis of excitation patterns, in the line of previous works on circular-core rotors and filaments (Keener, 1988;Wellner et al., 2002;Verschelde et al., 2007;Dierckx et al., 2012). The PDL concept provides not only a terminology for it, but also a way forward: solutions will need to be stitched together at the PDL interface, to elucidate their dynamics, much like was done before for wave front dynamics (Keener, 1986).
In the optical mapping experiments, the mechanical uncoupler blebbastatin was used, which however also affects the physiology of the cardiac cells. Therefore, the phase defect structures that we observed could be different from real cardiac tissue where no such uncoupler is present.
Furthermore, we used simple methods to find PDLs and WFs in the experimental data. As a result of representing the PDL as a spline curve, it is seen in Figure 13 that the wave fronts are not exactly touching the PDLs, while this is expected from theory. We will continue refining our numerical processing methods to make this image more consistent.
This study was also set in a fundamental science setting. Yet, it is inspired by the clinically relevant length and timescales: how do local dynamics self-organize in complex patterns with rotors and conduction block lines? In our opinion, contributions to the answer can be found from scrutinizing the patterns themselves, and building up a comprehensive understanding after the correct building blocks have been identified.
Finally, we suggest to try our methods to clinically relevant datasets, such as LAT maps derived from local electrograms. Existing algorithms could be tuned or extended, given that the classical PS was found in our simulations and experiments on a line of nearly discontinous phase.

Conclusion
In summary, the key ideas presented in this manuscript are: 1. Near conduction block lines, the phase of cardiac cells is close to discontinuous along a line, rather than along a point, as assumed in classical phase analysis; 2. Therefore, in excitable systems that have a forbidden zone in the inner region of their typical excursion (action potential) in state space, detected PSs are located at PDLs; 3. There is more than one definition of phase, and these different phases can be translated into each other depending on the goal and available data.
As a historical note, we find that the name of mathematician Bernhard Riemann is inscribed twice in our hearts: the heart is not only a Riemannian manifold due to anisotropy of wave propagation (Wellner et al., 2002;Verschelde et al., 2007;Young and Panfilov, 2010), but also features phase defect lines showing non-trivial Riemannian surfaces that organize the electrical patterns during arrhythmias.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by use of data (no new experiments).

AUTHOR CONTRIBUTIONS
HD and ET conceived the study. ET conducted optical mapping experiments. HD and LA ran the numerical simulations and performed experimental data analysis together with JQ, NN, and MV. All authors contributed to manuscript writing and internal reviewing.

FUNDING
ET was supported by National Science Foundation DCSD Grant 1662250. HD received mobility funding from the FWO-Flanders, Grant K145019N. LA was supported by a KU Leuven FLOF grant.