Directed Networks as a Novel Way to Describe and Analyze Cardiac Excitation: Directed Graph Mapping

Networks provide a powerful methodology with applications in a variety of biological, technological and social systems such as analysis of brain data, social networks, internet search engine algorithms, etc. To date, directed networks have not yet been applied to characterize the excitation of the human heart. In clinical practice, cardiac excitation is recorded by multiple discrete electrodes. During (normal) sinus rhythm or during cardiac arrhythmias, successive excitation connects neighboring electrodes, resulting in their own unique directed network. This in theory makes it a perfect fit for directed network analysis. In this study, we applied directed networks to the heart in order to describe and characterize cardiac arrhythmias. Proof-of-principle was established using in-silico and clinical data. We demonstrated that tools used in network theory analysis allow determination of the mechanism and location of certain cardiac arrhythmias. We show that the robustness of this approach can potentially exceed the existing state-of-the art methodology used in clinics. Furthermore, implementation of these techniques in daily practice can improve the accuracy and speed of cardiac arrhythmia analysis. It may also provide novel insights in arrhythmias that are still incompletely understood.


INTRODUCTION
One of the most effective ways to treat atrial and ventricular tachycardias is catheter ablation. In most of the cases ablation is guided by activation maps obtained from electroanatomical mapping systems. From these maps electrophysiologists need to precisely determine the mechanism of an arrhythmia-focal or reentrant-and assess the conduction pattern for a given patient to choose the proper ablation strategy. Performing it in complex substrates often confronts electrophysiologists with uncertainty (Delacretaz et al., 2001;Rostock et al., 2010;Kaiser et al., 2018;Martin et al., 2018). In these cases, the ablation procedure tends to be complex and time-consuming. This is particularly true for atrial tachycardias (AT) occurring after surgery or prior ablations (i.e., after ablation of persistent atrial fibrillation, Delacretaz et al., 2001;Deisenhofer et al., 2006;JaÏs et al., 2006;Patel et al., 2008;Rostock et al., 2010) and for scar-related ventricular tachycardias (VT) (Martin et al., 2018). If an incorrect target is ablated, not only will the patient not be cured, but new arrhythmias may be induced due to scarring (Chugh et al., 2005;Deisenhofer et al., 2006). In order to optimize catheter ablation, new tools for assessment of cardiac excitation patterns are needed to determine the underlying mechanism and to help identify the correct ablation target.
In the present study we propose a novel approach based on directed networks which allow the automatic determination of the type of cardiac arrhythmia (rotational or focal) and the characterization of important features of the excitation pattern which can be used for the automatic guiding of the ablation strategy. A network, in the most general sense, is a collection of nodes connected by links, which can represent diverse systems. Over the past 20 years, network theory has had many applications, ranging from biology to social sciences (Barabási, 2016). Examples include the PageRank algorithm (Brin and Page, 2012) for the World Wide Web which formed the basis of Google; determining the shortest route(s) between two places; modeling of molecules (e.g., fullerenes, Kroto et al., 1985), social networks (Borgatti et al., 2009), interactions of genes, proteins, metabolites and other cellular components (Barabasi and Oltvai, 2004;Barabási et al., 2011); the spread of diseases (Danon et al., 2011;Brockmann and Helbing, 2013); and many others. More recently, networks led to the development of novel diagnostic biomarkers in Alzheimer's disease, multiple sclerosis, traumatic brain injury and epilepsy (Stam, 2014). A network can be directed or undirected depending on if the connecting links have a direction from one node to the other. In spite of this variety, directed networks have not yet been applied to identify the sources of cardiac arrhythmias in the heart.
Directed networks naturally occur in the analysis of excitation patterns recorded by electrodes. When connecting discrete points of measurement in proximity to each other based on their local activation times (LAT), a directed network is created. This network of cardiac excitation appears suitable for directed network analysis. By applying network theory, conduction paths can be identified in a new and different way based on local activation times and by taking the physiological conduction velocity into account. During the analysis, the algorithm identifies potential ablation targets such as rotational activity, spreading from electrode to electrode creating a closed loop, or focal activity, manifesting as a divergence of excitation from a given point (region). We refer to this method as directed graph mapping (DG mapping). In graph theory, very efficient methods have been developed to find closed loops in directed networks. Using these methods, one can easily find all possible loops in these data within mere seconds. Since this approach analyzes all possible loops automatically, it forms a robust method even in the presence of noise or incorrect electrode recordings, making it much more reliable than current existing methods. Furthermore, it allows the determination of additional properties of excitation as well, which can be essential for the characterization of the arrhythmia. By using directed networks and DG mapping, we believe that a more reliable, faster and fully automatic analysis of activation patterns can be performed with a higher accuracy than in current daily practice.
The goal of this study is to demonstrate the wide applicability of directed networks to the heart for each driving mechanism of cardiac arrhythmias in both the atria and the ventricles. Therefore, we tested the accuracy of DG mapping in in-silico (ventricular) models of functional and anatomical reentry and focal activity. To determine the accuracy of DG mapping in the atria, we analyzed 31 clinical cases of atrial tachycardia. Regular AT is a clinical tachyarrhythmia in which the operator can be sure about the location of the tachycardia since ablation of the correct target almost always results in immediate success. Therefore, AT was used as the gold standard for validating DG mapping in a clinical setting. In addition, DG mapping was compared to phase mapping (Gray et al., 1998) via in-silico simulations, a widely used technique for detecting the center of a rotor.

MATERIALS AND METHODS
In the next sections the general protocol of DG mapping is explained according to the flowchart given in Figure 1.  Figures 2, 3 demonstrate DG mapping on a simulated and a clinical example.

Input
The DG-protocol can be applied to a wide range of in silico, experimental and clinical models of arrhythmia with different types of electrode systems (e.g., basket electrode system, intramural needle electrodes, high density grid data, unstructured electrodes systems, etc.). In the current study, this will be demonstrated with in-silico generated and clinical data.

In-silico Generated Datasets
All simulations were performed using the TNNP-computer model for human ventricular cells (ten Tusscher and Panfilov, 2006) utilizing the explicit-Euler integration scheme (Vandersickel et al., 2014) on a GeForce GTX 680 and a GeForce GTX Titan, with single precision. The following different scenarios were simulated: (1) Functional reentry was simulated in 2D (in a domain of 512 by 512 grid points with interspacing of 0.25 mm) and 3D (a simplified model of the human ventricle and an anatomically accurate model of the human ventricle Tusscher et al., 2007). (2) Anatomical reentry was also simulated in 2D and 3D (anatomically accurate model of the human ventricle). In both scenarios, the S1S2-protocol was applied to obtain rotational activity (Tusscher and Panfilov, 2003). (3) Focal activity was simulated in 2D and 3D (anatomical model of the human ventricle) by applying 3 stimuli at 3 different locations of 500 ms each. All simulations were performed for a duration of 20s. In all simulations, the rotors were stable in space and time.
For each different setup, we implemented either 64 surface electrodes (mimicking 64 electrode-basket catheters, Narayan et al., 2012), 256 surface electrodes with an interspacing of 0.8 mm (mimicking experimental grid sizes, de Groot et al., 2010) or 500 intramural electrodes (in the 3D anatomical model) in analogy with the experimental setup by Taccardi et al. (2005) In Figure 2, an example of a rotor with 64 electrodes is shown. For these electrodes, we computed local unipolar electrogram as follows: FIGURE 1 | Illustration of the work flow of the DG mapping tool. As input, the DG-tool requires for a given setup of electrodes the LAT values with the corresponding XYZ coordinates, which can be extracted from either simulation studies or a clinical setup. The input is then processed as follows in the DG-tool as presented in the flowchart. Next, we apply a loop-finding algorithm to detect cycles in the network. If cycles are not detected, we locate the source of focal activity. In case cycles were detected, the loops are merged and its center is determined. At the end, the output is visualized. In case of a focal source, arrows pointing away from a (group of) node(s) are shown, while for reentry, arrows will be plotted to visualize the reentry path.
where t is time, x is the distance to the tissue, V is the transmembrane potential and r is the integration parameter over the tissue. The XYZ-coordinates of the selected electrodes were also stored for further analysis. The LAT of each electrode was determined by taking the steepest negative slope (−dV/dt) of the calculated unipolar electrogram, see also Figure 2A.
This coincides with the upstroke of the action potential (i.e., the true moment of activation) (Spach et al., 1979;Spach and Dolber, 1986).

Clinical Datasets
Between April and August 2017, 29 patients undergoing ablation of symptomatic ATs at AZ Sint-Jan Bruges Hospital were enrolled in the study, resulting in 31 activation maps (30 left atrium, 1 right atrium). The study was approved by the local ethics committee of AZ Sint-Jan Hospital Bruges. High density (> 300 points) endocardial mapping of ATs was performed using a single-electrode mapping and ablation catheter with a distal 3.5 mm tip and three 1 mm ring electrodes (THERMOCOOL SMARTTOUCH Biosense-Webster Inc., Diamond Bar, CA, USA). These high density maps covered the full atrium. The bipole of a decapolar coronary sinus (CS) catheter was selected as reference for activation mapping (i.e., peak of CS = 0 ms). The following settings for activation mapping were applied: mapping window set to tachycardia cycle length minus 10 ms and centered at the 0 ms reference. Usually the activation map window is set in this way with the aim to cover the entire tachycardia cycle length during mapping (Del Carpio Munoz et al., 2010). This mapping window is a filter criterion used during continuous mapping that compares LATs between two consecutive beats, but only if the difference in LAT does not exceed 10 ms are the data then acquired. This filter enables correct and accurate data acquisition in order to make a consistent activation map. The other settings were minimum contact force of 4 g, LAT stability of 10 ms, respiratory gating enabled and color fill calibrated at 5. Bipolar scar threshold was defined at 0.05 mV (Anter et al., 2015), and EGMs with bipolar voltages lower than this cutoff were therefore automatically tagged as scarring (gray zones) on the activation maps. Automated and continuous acquisition of points was performed by the CONFIDENSE mapping module (Carto 3 v. 4, Biosense Webster Inc.) using the novel hybrid LAT annotation method (LATHybrid) (Pooter et al., 2018). Each AT case was analyzed offline by DG mapping after exporting all local  Frontiers in Physiology | www.frontiersin.org activation times (LATs) and the corresponding 3D coordinates. In Figure 3A, an example of the left atrium is shown, with the corresponding LAT map and annotated points.
The tachycardia mechanism was confirmed when ablation resulted in sinus rhythm or in conversion of a second tachycardia. In case of multiple hypotheses of the AT mechanism, the hypothesis which agreed with the ablation endpoint was chosen.

Directed Graph Mapping Protocol
This section explains the DG mapping algorithms, as shown in the blue panels in Figure 1.

Determine the Neighbors in a Given System
First, for a given configuration of electrodes, possible neighbors for each electrode are determined. These neighbors cover all possible paths where the wave can travel to, starting from a certain electrode. For regular grids, the neighbors are found by setting a spheric distance around a single point. Hence, a single point incorporates up to 8 neighbors in case of the 2D grid (see Figure 2B) and up to 26 neighbors in case of a regular 3D grid. For an irregular configuration of electrodes, like the clinical AT cases, Delanauy triangulation is applied to determine for each electrode its possible neighbors (see Figure 3B).

Creating Network of Cardiac Excitation
We chose a certain time t. Starting from this time, we find LAT 1 , ...LAT n which are the first LAT larger than t for each electrode in our system of n electrodes. We then draw arrows as follows. Suppose electrodes 1 and 2 form a pair of neighbors. Assume electrode 1 has LAT 1 and electrode 2 has LAT 2 , with LAT 2 > LAT 1 , meaning the difference between the two electrodes is δLAT = LAT 2 − LAT 1 > 0. We allowed a directed link from electrode 1 to 2 if : In this equation, CV min , CV max , and d represent minimal conduction velocity, maximal conduction velocity and the euclidean distance between the two electrodes, respectively. For the simulated examples for ventricular tissue (2D and 3D) we took CV min = 0.2 mm/ms, and CV max = 2.00 mm/ms. For the clinical AT cases, CV min was set at 0.08 mm/ms, according to the lowest physiological conduction velocity in human atria determined by Konings et al. (1994), CV max was set to maximal 2.0 mm/ms (Harrild and Henriquez, 2000). In Figures 2C, 3C, the directed arrows from a single electrode are shown. Once this first graph was created, a second graph at a time t + δt was created in exactly the same way as the first graph. We set δt = 40 ms. Finally, these two graphs were merged, whereby arrows of the second network were added to the first network if the LAT of the node where the arrow originates from was the same. This was necessary as in the first network, no closed cycles will be present, which represent the rotational activity of the arrhythmia, and they are exactly the arrows of the second graph, which will create cycles in the network. The resulting graph is the final directed network. δt = 40 ms is an arbitrarily chosen value to make sure the network indeed forms closed loops, but the algorithm can work as well for other values of 0 < δt < CL/2. However, one cannot make δt too small, as otherwise the first and second graph might be equal. For example, in Figures 2D, 3D, the complete network is shown for a simulated case and a regular AT.

Rotational Activity
Once the network is created, any type of rotational activity can be found by detecting cycles in the network. A cycle is a closed directed walk without repetition of nodes. In order to find the cycles, a standard breadth-first search algorithm was used. Since the constructed network generally turns out to be rather small and very sparse, this can be done very efficiently. It turns out that detecting all (smallest) cycles through each node can be done almost instantaneously. We ran theoretical simulations on networks with 1,000,000 nodes, and even in these cases all cycles were found in the range of seconds. Clearly, the physical bounds on the number of electrodes that can be placed will be more limiting than the computational work that is needed to process the data. In Figures 2E, 3E, the resulting cycles of the network of a simulated rotor and a regular AT case are shown.
In order to find the core of any type of rotational activity, we looked for the smallest cycles in the network and computed the geometric center. This was performed by grouping all found cycles based on their proximity to the geometric center. If the centers lie closer to each other than a specified threshold, the cycles were considered to belong to the same core. In this study, we took 1 cm as threshold, as we estimated that the cores of the reentry loops considered in this work were always apart > 1 cm. Afterwards there was an optional pass which merges bundles of cycles if they shared nodes. Finally, the centers of each bundle were defined as the core of rotational activity. In Figures 2F, 3F, the selected cycles are shown.

Focal Activity
Focal activity was detected as a source node, i.e., a node which has a non-zero out-degree, and an in-degree equal to 0. These can be found immediately by doing a single pass over all nodes. Then, the LATs were bundled in certain intervals to reduce the inter-variability in the LAT values. Afterwards, we reconstructed the network with these bundled values. We then checked if regions with only outgoing arrows were present. The middle of these regions corresponds to the source of the focal activity. In Figure S1, we have repeated all the previous steps for a simulated focal source.

Additional Features of Network Theory
In addition to finding rotational and focal activity, we derived additional properties of the network.

Region of Influence
For each network containing reentrant circuits, we can determine a "region of cycles" and a "region of influence." The region of cycles contains all nodes (electrodes) which are part of cycles for a particular reentrant circuit. Second, for each non-marked point we can determine the closest "region of cycles" in terms of network arrival time distance and relate it to that region. As a result, for each point we can determine which source excited it. This is called the "region of influence." In order to construct the region of influence, the following algorithm was implemented. For a given network, all n cores were determined, c 1 , ..., c n . For each core, we first determined all nodes which are part of cycles of the network (C 1 , ..., C N ), i.e., the regions of cycles. Then, each node was added to the core c i to which it had the shortest path to one of its nodes in C i . In this way, each core is assigned a region of influence.

Wave Averaging
Another application of the constructed network is wave averaging to interpret the cardiac excitation pattern. In general, the outgoing arrows of each node were averaged, and only this average arrow was kept for the visualization. In more detail, the following steps were taken in the wave averaging algorithm. First, each LAT-node was projected on the geometry (mesh) of the atrium. Second, each arrow of the directed network was projected by dividing the arrow in 4 equal parts and projecting these parts on the geometry. The begin and endpoints of these arrows form new nodes which were added to the existing nodes. Then, for each node n on the geometry, each directed arrow starting from this node as well as each connection of each node within 1 cm from the original node n was averaged. The collection of these averages was then plotted on the geometry.

Phase Mapping Protocol
LAT values were used to construct the excitation patterns in phase-space. First, a sawtooth wave, with amplitude ranging from −π to π is constructed based on these LAT values. Afterwards, values are adapted with their 2π equivalent within the range of −π to π in phase-space. Next, in both x and y directions, the phases were derived and a linear combination with the Sobel (we also tested the Nabla) kernels to detect the singularities was applied. This protocol was previously presented (Bray et al., 2001;Bray and Wikswo, 2002a,b;Umapathy et al., 2010;Tomii et al., 2016). However, based on the properties of the ECGs of the simulation, we made use of this sawtooth wave instead of the more regular Hilbert transform of the ECG signal, as this is does not make any difference for regular signals-see Kuklik et al. (2015). In 3D, the heart was sliced in 3 orthogonal directions and the protocol was applied on each slice. However, as the shape of the ventricle model is complex, only grid points with complete circumference in the heart were taken into account, so convolution did not result in false positives on the edges (Van Nieuwenhuyse et al., 2017). However, this did not result in detection of the filament of the rotor as the density (500 intramural points) was too sparse. We therefore calculated the phase singularities on the surface of the tissue and detected eventually the phase singularities of the spiral in 3D. A binary detection threshold was applied to the convolution (Tomii et al., 2016), set to 95% of the maximal detected value in phase-space.

Introducing LAT Variation
In the clinical setup, identification of LAT either by automated algorithms or manual annotation by operators can vary due to several factors such as accuracy of the detection algorithms, operator experience, signal quality and noise (El Haddad et al., 2013). Therefore, we included LAT variation in our analysis, and compared the accuracy of DG mapping with phase mapping. In order to obtain LAT variation in the simulations of functional reentry, random Gaussian noise was added with standard deviations σ = 5, 10, 15, 20, 25, 30 ms in the simulation of functional reentry with a configuration of 64 and 256 electrodes. We divided the activation times obtained during a simulation in 25 different frames with 520 ms separation to exclude any overlap in activation times. For each frame, we randomly added Gaussian noise 1,000 times, so in total, we compared 25,000 different frames per LAT variation σ . The center of the rotor was detected through DG mapping and phase mapping. For DG mapping, the geometric center of all cycles belonging to the same core was computed. Afterwards, the median value was taken as the true center of the rotor. In addition, only the center with the highest number of cycles was taken into account. We classified the outcome as correct if only one single core was found within 1 cm of the true core. The incorrect diagnosis was classified in 3 different types: incorrect cores (i.e., cores outside a 1 cm radius of the true core) in combination with the correct core (error type 1), only incorrect cores (error type 2), or no cores (error type 3). For the percentage correct diagnosis p, we computed a 95% confidence interval via p ± 1.96SE where SE is obtained from a robust sandwich estimator (Thomas Lumley, 2015; R Core Team, 2017) that accounts for the correlation structure (i.e., the 1,000 replicates within one time frame are expected to be correlated). In the Supplementary Material, we also simulated noise from the skewed lognormal distribution to study the robustness of the methods for different types of noise distributions. In addition, we also presented the outcome as a function of the distance from the true core (instead of taking 1.0 cm).

In-silico Models of Functional and Anatomical Reentry and Focal Activity
The accuracy of DG mapping was tested in different in-silico models as described in the methods section. First, for functional reentry (see Figure 4A), we simulated a 2D rotor with a configuration of 64 electrodes (A1) and 256 electrodes (A2). In 3D, functional reentry was induced in a simplified model of the ventricles with 64 surface electrodes (A3) and in an anatomical model of the ventricles with 500 intramural electrodes (A4). In all four setups, DG mapping was able to accurately detect functional reentry and correctly determine the location of the core of the rotor for the entire length of the simulation (20 s duration). The smallest cycle and corresponding core are shown in yellow for each setup. Second, DG mapping was validated in two models of anatomical reentry ( Figure 4B): a 2D anatomical circuit with 64 electrodes (B1) and a 3D anatomical reentry with 500 intramural points in the model of the ventricles (B2). In both models, DG mapping correctly identified the reentrant path around the obstacles for the entire length of the simulation (20s). The shortest reentry loops are again depicted in yellow. Third, focal activity was simulated in 2D (64 electrodes) and 3D (500

Clinical Dataset
To establish proof of concept in the clinical setting, we retrospectively and blindly analyzed 31 cases of regular atrial tachycardia (AT). For clarity, in Figure 3, all the steps of the DG mapping protocol were demonstrated on an AT case of a localized reentry. In general, the atria have a complex structure. In case of reentry during AT, the electrical waves circle around obstacles such as the valves, the veins or scar tissue, creating a (sustained) reentry loop. Ablation aims to terminate the reentry loop so that the circular electrical conduction can no longer be sustained.Therefore, it is important to precisely determine the location of the activation pathway. The accuracy of DG mapping was compared to the standard diagnosis, i.e., type of arrhythmia and location of the circuit/focal activity as determined by the electrophysiologist (EP) based on the activation map and the ablation result. The overall results are summarized in Figure 5. Out of 31 cases, 20 were due to macro-reentry, 6 due to localized reentry and 5 due to focal activity-see also the Table S1. In 9 cases with reentry, the operator was not sure about the reentry mechanism purely based on the LAT activation map, formulating several hypotheses. The gold standard was taken as the diagnosis FIGURE 5 | Accuracy of DG mapping in clinical AT. The gold standard was compatible with focal source (5), localized reentry (6), and macro-reentry (20). Macro reentry was categorized in reentry around the right veins (RV), the left veins (LV), mitral valve (MV), around RV + LV, around RV + MV, around LV + MV, and other types of reentry (e.g., in the right atrium). In brackets the accuracy of DG mapping is given.
matching the ablation endpoint. Compared to this gold standard diagnosis, DG mapping identified the exact same mechanism and location in 28 out of 31 cases (90.3%, 95% exact binomial confidence interval 74.2% -98%). In 3 out of 31 cases, the diagnosis of DG mapping did not fully match with the gold standard. In 2 cases of double loop reentry (cases 6 and 14), DG mapping identified only one single loop. In the other case (case 22), the mapping data indicated focal tachycardia, whereas DG mapping identified localized reentry at the same location. However, in all 3 cases, DG mapping would have pointed to the correct ablation target, meaning that DG mapping correctly identified the ablation target in 31/31 cases.
Representative cases are shown in Figure 6. Panel A depicts a macro-reentrant AT around the right pulmonary veins in the LA conducting over the roof. Ablation of the roof resulted in prompt termination of the AT. Blinded analysis by DG mapping revealed a selected loop at the same location (middle panel). Panel B shows a localized reentry at the anterior wall, rotating around local scar tissue. Ablation from the scar to the mitral valve terminated AT. DG mapping (middle) as well as wave averaging (bottom) identified the same location of the localized reentry. In panel C, activation mapping and ablation conformed with focal tachycardia at the septum. DG mapping (in the absence of loops) pointed to focal activity as well (middle panel). We also tested the wave algorithm for each case, as shown in the bottom panels of Figure 6. Wave averaging was for all the cases compatible with the results of the DG mapping. Representative examples are shown in Figure 6: macro reentry (A), localized reentry (B), and focal activation (C).

Comparison With Phase Mapping Under LAT-Variation
In the clinical setup, LATs can vary due to several factors such as accuracy of the detection algorithms, operator experience, signal quality and noise (El Haddad et al., 2013). Therefore, the performance of DG mapping was compared to phase mapping in the model of a single rotor with 256 electrodes, now by adding Gaussian white noise to the LATs (Figure 7, upper panels). Overall, we observed that DG mapping retains its accuracy to detect rotors at increasing noise levels, whereas phase mapping becomes less precise (middle panel of Figure 7): for small variation levels (5 ms), DG mapping is 100% accurate, while the accuracy of phase mapping decreased to 74.17%. For 15 ms, phase mapping became highly unreliable (accuracy of 30.22%) while DG mapping had an accuracy of 95.49%. For 20 ms, this difference was even more pronounced: DG mapping maintained an accuracy of 81.19% while the accuracy of phase mapping dropped to 1.08% (p-value < 0.0001). Moreover, in case of incorrect diagnosis (lower panels), phase mapping detected extra false cores (type 1 error) whereas in the DG method incorrect diagnosis was due to no detection of the core (type 3 error). Noise analysis was repeated for the 2D model with 64 electrodes. A similar trend was found with DG mapping being more accurate (91%, 83%, 68%, for a noise level of 10, 15 and 20 ms) vs. phase mapping (75%, 51%, 21%, respectively). All of these differences were highly significant (p-value < 0.0001), see also Table S2. We also varied the distance to the true core for which we retained a diagnosis as correct, as shown in Figure S2 (and see also Figure S3 for more explanation). As explained in the supplementary material, due to the discrete nature of phase mapping, we can only compare phase mapping and DG mapping above a certain threshold. In these cases, DG mapping always exceededs phase mapping. We also tested the effect of the underlying distribution of the LAT values. Modeling LAT variation with the skewed lognormal distribution did not alter the conclusions of the results-see also Figure S1. Finally, to test the specificity of DG mapping, we applied DG mapping in a point stimulation model without functional reentry (256 electrodes, LAT variation ranging from 0 ms to 30 ms). In these cases, DG mapping never identified any rotors, resulting in a specificity of 100%.

Region of Influence
Describing cardiac excitation as a network allows the extraction of additional information. Besides wave averaging, DG mapping allows the identification of the spatial region which is excited by a certain source. In case of normal excitation, a single source (sinus node) excites the whole medium. However, in case of an arrhythmia with multiple sources, each source excites a given region, which we call the region of influence. We hypothesized that DG mapping, by containing complete spatiotemporal information, can determine this area of influence. This concept was evaluated in 2 different setups (see Figure 8A). We determined the region which contained the electrodes belonging to cycles (region of cycles) as well as the region of influence. Obviously, for a single rotor, the region of influence spans the entire set of electrodes. For 4 different rotors, one can observe an area of influence for each given source.

Main Findings
In this paper, we demonstrated that a directed network can be used to describe electrical activity in the heart in order to find the mechanism of the arrhythmia. This novel method is robust, fast, general, accurate and can be applied to a wide range of in silico, experimental and clinical models of arrhythmias. First, we showed that DG mapping can find functional reentry, anatomical reentry and focal activity in in-silico ventricular models of the heart (see Figure 4). We tested using intramural electrodes, 64-basket electrodes and regular grids with different numbers of electrodes (64-256). Second, we tested DG mapping on 31 clinical cases of regular AT-see Figure 6. Compared to the gold standard, DG mapping identified the exact same mechanism and location in 28 out of 31 cases, whereas it identified the correct ablation target in all 31 cases. These results suggest that DG mapping could potentially lead to improved treatment of tachyarrhythmias based on stable sources.

Network Theory
To our knowledge, so far only limited research has focused on network theory to understand cardiac arrhythmias. In the study by Zahid et al. (2016), undirected networks were used to find the minimal number of nodes which need to be ablated to separate two regions in the heart. This region was then proposed as the optimal ablation site. In the study by Tao et al. (2017), the authors showed that ablation of persistent AF is associated with improvement in both local and global connectivity within the communication networks. However, in neither of the above studies, excitation was interpreted as a directed network. Zeemering et al. (2013) applied a directed network to describe AF by accumulating multiple time frames. However, in contrast to DG mapping, this methodology precluded the possibility of detecting rotational activity and does not represent the actual wave excitation. Also in Richter et al. (2012), AF was described as a directed network via the use of sparse modeling for the estimation of propagation patterns in intracardiac atrial fibrillation. However, it is not clear how rotational activity can be detected from the obtained networks, and it would be of interest to uncover the cycles. Similarly, in Alcaine et al. (2016) and Luengo et al. (2018), directed arrows are created based on the concept of Granger causality between different signals instead of the LATs. This could form an alternative way to create the directed network if deriving the LATs from the signals is not feasible. In a model of chronic atrial ventricular block, we used directed networks for the first time to determine the mechanism underlying Torsade de Pointes . However, again, we did not use it to fully describe the electrical excitation as in this work. Therefore, to our knowledge, this is the first study where directed networks were used to FIGURE 7 | Performance of DG mapping and phase mapping (PM) while adding Gaussian noise with an increasing standard deviation (ranging from 0 to 30 ms) in a single 2D rotor model with 256 electrodes. The upper panels show representative activation maps for different levels of noise. The middle panel shows the performance of DG mapping vs. phase mapping for these different noise levels. The bottom panel shows the type of errors in case of failure for both methods. Error type 1 is a detection of a false core in addition to the correct core, error type 2 is only a detection of a false core and error type 3 is no detection of cores at all. describe electrical excitation to extract the mechanism of the arrhythmia, building on our previous work in the CAVB dog .

Advantages of DG Mapping
First, we showed that DG mapping could be used to reliably detect rotational activity even after adding LAT variation. For instance, in the model presented in Figure 7, we found that in case of 15 ms standard deviation of noise, phase mapping was only 30% accurate, while DG mapping was still 96% accurate. This difference in accuracy can be explained by the holistic nature of DG mapping. In contrast with phase mapping, in the presence of a number of electrodes with a wrong LAT annotation, DG mapping can still identify the correct location of the rotor based upon the other electrodes. In addition, DG mapping also takes into account the number of cycles which are found for each rotor (with a higher number of loops indicating a higher likelihood of an actual rotation). In contrast, phase mapping finds phase singularities locally. Therefore, a misplaced LAT can easily give rise to false positives/negatives, potentially resulting in incorrect clinical decisions. Also, DG mapping takes into account the conduction velocities of the excitation, which corrects for non-regular spacing of the electrodes, a feature which phase mapping lacks.
Second, DG mapping automatically detects focal activity, a feature which phase mapping lacks. Also, depending on the size of the obstacle and the spacing of the electrodes, phase mapping cannot detect anatomical reentry. In addition, phase mapping might indicate anatomical reentry around a small obstacle as functional reentry. Very recently, another methodology was described by Oesterlein et al. which can also automatically detect anatomical reentry (Oesterlein et al., 2018). The approach uses a different methodology based on integral measures: determination of activated area and its relation to the cycle length of the arrhythmia, while our DG mapping directly analyzes the local propagation of the excitation wave. It would be interesting to compare this method with DG mapping, especially in clinical settings. Also, other techniques exist to help uncover the mechanism for atrial arrhythmias like the ripple map (Linton et al., 2009) (AT) or the retro mapping technique (Mann et al., 2019) (AF). However, both methods still require a manual interpretation of the novel types of maps which were obtained by these techniques. Also very recently, the STAR technique was developed (Honarbakhsh et al., 2019). The STAR methodology aims to identify the predominant wavefront direction by displaying it on the STAR maps. Although it does not focus on identifying a particular arrhythmia mechanism, it helps the operator to determine the mechanism from the STAR maps.
Third, DG mapping offers additional features which can be derived from the directed network. DG mapping can determine all electrodes belonging to any cycle which are part of the same rotational activity (see Figure 8A) and detect for each rotational core its region of influence ( Figure 8B). This offers the possibility to detect all electrodes activated by a specific rotational activity and could detect the dominant driving source of the arrhythmia as a primary target of ablation. Also, the wave averaging technique allows the creation of maps of the wave propagation, which can provide additional guidance during catheter ablation.
Finally, another advantage is that DG mapping is universal. It can be applied to any type of recording system from which LATs can be derived, with varying number of electrodes, the inter-electrode distance or site of recording, as shown in

Clinical Implications
As shown in this paper, DG mapping can be of added value in the ablation of regular AT. Despite improvements in activation mapping (RHYTHMIA by Boston Scientific, Coherent mapping system by Biosense Webster), interpretation of activation maps remains challenging and operator dependent (Gerstenfeld and Marchlinski, 2007;Kaiser et al., 2018). We demonstrated that DG mapping automatically identified the same mechanism as the electrophysiologist (EP) in 28/31 cases of regular AT, but did find the correct ablation target in 31/31 cases. Moreover, in 9 cases with reentry, the operator was not sure about the mechanism based on the LAT activation map and formulated several hypotheses (see Table S1). Therefore, DG mapping could aid physicians in finding the correct diagnosis according to the ablation target. Currently, in case of doubt, the operator can perform entrainment mapping whereby the post pacing interval (PPI) is compared to the tachycardia cycle length (TCL) to localize or confirm the correct reentry circuit (Knight et al., 1999). Furthermore, DG mapping also automatically detects focal activity and its location, making it a complete diagnostic tool for AT. Compared to the standard assessments, DG mapping is robust, fast and operator independent. Therefore, DG mapping would remove the manual interpretation of the (experienced) operator. The wave averaging technique allows the creation of maps of the wave propagation, which can provide additional guidance during catheter ablation. Another advantage is that DG mapping is instantaneous and could therefore shorten the ablation procedure.
Another important potential application could be AF. AF is often referred to as the most common arrhythmia in clinical practice, with an estimated prevalence of 2%, and is associated with a fivefold and twofold higher risk of stroke and death, respectively (Zoni-Berisso et al., 2014). Catheter ablation of AF yields moderate success rates (Brooks et al., 2010;Weerasooriya et al., 2011;Weiss et al., 2016), which is related to the lack of understanding of AF mechanisms. Different mechanisms for AF have been described such as focal activation, dissociated activity or stable rotors (Allessie and de Groot, 2014a,b;Narayan and Jalife, 2014a,b;Chen et al., 2015;Kirchhof et al., 2016). Currently, both researchers and electrophysiologists rely on activation mapping or phase mapping for the analysis of AF. Recently, initial studies suggested good outcomes after ablation of rotors guided by phase mapping (Narayan et al., 2012(Narayan et al., , 2014. However, new studies have emerged contradicting this study (Buch et al., 2016;Gianni et al., 2016;Mohanty et al., 2016). It was shown that phase mapping easily generates false positives (Vijayakumar et al., 2016;Kuklik et al., 2017), especially due to LAT variations on the signals and large inter-electrode distances. In our study, phase mapping showed similar results when adding LAT variations, whereas DG mapping maintained a high accuracy. It remains to be seen whether DG mapping will offer new insight in AF mechanisms; however, the holistic nature of the method (as explained in the advantages of DG mapping) might overcome the problems with phase mapping, as currently used in the clinic.
In conclusion, translating cardiac arrhythmias into directed networks as described in the current work offers the beginning of a new area of research. There exists a whole range of different algorithms in network theory (e.g., edge density, centrality measures, etc. White and Borgatti, 1994;Newman, 2010;Holme, 2015;Sizemore et al., 2018), which can possibly be applied to the constructed networks to increase our understanding of cardiac arrhythmias.

Limitations
As this paper was a proof-of-concept, many different settings are not yet tested. For example, it remains to be tested how DG mapping will characterize cardiac excitation in more complicated settings with multiple meandering rotors, including wavebreaks, or in complex fibrotic tissue. A limitation of DG mapping is that it requires at least one full cycle (as DG mapping can only find closed loops) of a circular rotation, while phase mapping finds phase singularities instantaneously. A limitation of DG mapping is that it requires a degree of stability of the cardiac excitation pattern since a full cycle of activation is required for DG mapping. Therefore it remains to be seen in the clinical setting if DG mapping can advance the understanding of more complex arrhythmia such as AF, VT and VF. For these cases, DG mapping requires the arrhythmia to be mapped first, which is not always possible (e.g., due to hemodynamic instability in VT/VF etc.). Future studies are needed to further evaluate DG mapping in different types of arrhythmias.

DATA AVAILABILITY
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by ethics committee of AZ Sint-Jan Hospital Bruges. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
NV coordinated the project, set up the collaboration, performed the computer simulations, and developed the DG-mapping protocol. EV developed and optimized software for the study of DG-mapping in simulations, all clinical AT cases, and analyzed simulations with phase mapping. JG and NVC developed the software to analyze the networks created by NV and EV. ME provided all data and analyzed the clinical data. JD performed the statistical analyses. TS contributed in analyzing the clinical cases. AD contributed in writing the paper. MD analyzed all 31 clinical cases, performed all ablations which were analyzed in this work, and significantly contributed to all parts of the paper. AP contributed to all parts of the paper and co-developed a part of the protocol (region of influence). All authors contributed and collaborated in writing and improving all aspects of the paper.

ACKNOWLEDGMENTS
This manuscript has been released as a Pre-Print at Vandersickel et al. (2019).