Phase Gradient Divergence for the Quantitative Detection of Focal Activation Events During Cardiac Fibrillation

Atrial fibrillation is the most common arrhythmia. Spiral wave and focal activation (FA) are known to play an important mechanistic role in the generation, sustenance, and termination of tachyarrhythmia. However, to date no quantitative method of detecting FA under complex excitations has yet been established. In this study, we propose such a detection method of FA by calculating the divergence of the spatial gradient of a phase map, which identifies the phase of the excitation cycle at each point of the heart as derived from measurement signals. Next, to verify the accuracy of the proposed method, we conducted a membrane potential measurement experiment using an excised porcine heart (n = 1). By comparing the conventional and proposed methods for 126 instances of FA, we found that the proposed method showed improved detection accuracy, with precision, sensitivity, and F-measure values of 0.45, 0.84, and 0.58, while conventional method showed 0.04, 0.26, and 0.08, respectively. The proposed method, which uses the divergence of the phase gradient to predict FA, shows potential to suppress false positives that are observed in the conventional method, and it can more accurately detect FA than conventional methods.


Introduction
Tachyarrhythmia is a condition in which the normal heartbeat is disrupted and the heart beats faster than normal cycle due to abnormal excitations, resulting in atrial fibrillation (AF) and/or ventricular fibrillation (VF). AF increases the risk of serious and fatal complications, such as stroke and heart failure [1,2], while VF is a cause of sudden cardiac death [3]. It is therefore necessary to establish an effective treatment for these tachyarrhythmias.
Basic research has been conducted to clarify the mechanisms underlying tachyarrhythmia with the aim of establishing an effective treatment. Recently, it has been suggested that spiral waves (SWs) and focal activation (FA) play a crucial role in the generation, sustenance, and termination of tachyarrhythmias [4,5]. Therefore, to control tachyarrhythmia, it is necessary to accurately understand these complex excitations that occur during arrhythmia. Moreover, to achieve this, it is essential to develop stable detection methods for SWs and FA.
Various attempts have been conducted to develop SW detection methods [6,7]. For example, a phase variance analysis proposed by Tomii et al. may be able to suppress false positives that are frequently observed in conventional methods used to detect SWs. Thus, this method may permit enabling stable detection of SWs [8].
Similarly, there have also been attempts to develop a reliable detection method for FA [9,10]. However, these detection methods have not been sufficiently evaluated because most of these methods use clinical data in which the true presence or absence of FA is unclear. On the other hand, other studies have quantitatively evaluated methods to detect FA using computer simulations [11][12][13]. In these studies, an activation map is first derived by spatially mapping the time until the next activation based on the detected activations measured using electrode signals. Subsequently, FA is quantitatively detected using numerical methods, such as spatial gradient and divergence analysis, or the directed graph analysis of the activation map. However, it is known that the 1 3 conduction velocity (CV) vector field calculated using the activation map is discontinuous near the wavefront (WF) [14]. It is therefore expected that numerous false positives will result from applying this method to complex and aperiodic excitation during cardiac fibrillation, where multiple WFs occur simultaneously. Therefore, to date, no quantitative detection method of FA under arrhythmia conditions has yet been established.
In this study, we propose a quantitative detection method of FA under complex and aperiodic excitation. The proposed method involves calculating the divergence of the spatial gradient of the phase map, which spatially maps periodic excitations as circular variables. We verified the accuracy of the proposed method by comparing with a conventional method by conducting ex vivo experiment using porcine atria. Figure 1a shows the problem in the conventional detection method and an overview of proposed method. In the conventional method, an activation map is plotted by mapping the time until the next activation at each point measured by electrode measurement signals. Subsequently, FA can be quantitatively detected using spatial gradient and divergence analyses or the directed graph analysis of the activation map. However, due to discontinuous values around the WF in the activation map, the direction of CV vectors is reversed near the WF (Fig. 1a, upper panel). This is expected to cause high divergence values near the WF even though FA does not occur; this can therefore be a cause of false positives. In this study, we propose a detection method of FA using a phase map instead of an activation map. A phase map is a spatial mapping of the phase assigned to the periodic excitation of myocardial cells, and the phase value at each point is a circular variable ranging − π to π. It is expected that the discontinuous boundary near the WF that occurs when using the conventional method will be eliminated using a circular variable for gradient calculations, and this in turn will cause false positives to be suppressed (Fig. 1a, lower panel). Figure 1b shows the analysis flow of the proposed method. Here, after calculating and normalizing the spatial gradient of the phase map derived from the cardiac excitation, the divergence is calculated as follows: The analysis flow of the proposed method. After deriving the phase map from heart excitation measurement, we calculate and normalize the spatial gradient and then compute the divergence. Additionally, an excitation wavefront (EWF) map is multiplied by the divergence map so that only the values at the FA occurrence time are enhanced. Finally, we binarize the resulting map and labeled the area used to detect FA where P is the normalized phase gradient. In addition to calculating the divergence, an excitation wavefront (EWF) map is derived by transforming the phase value of each pixel of the phase map as follows:

Proposed Method
where θ is the phase value, θ act is the phase value corresponding to the WF, and σ is the parameter governing the spread of the exponential function. As the phase gradient continues to spread concentrically around the FA position as long as it occurs, we multiply the divergence of the phase gradient by this EWF map so that the values are higher only during the time. Following this calculation, we binarized this map using a fixed threshold and labeled the positive area. Finally, we calculated the center of gravity of each labeled area in the spatial and temporal directions and determined the position and time of FA occurrence.

Experiment
We conducted the following experiment to compare and evaluate the accuracy of FA detection using the conventional and proposed methods. Figure 2 shows an overview of the experimental system.

Porcine Model of AF Under Acute Intra-atrial Pressure
The protocol followed in this study was approved by the Institutional Animal Care and Use Committee of the University of Tokyo. One porcine was anesthetized and its heart was quickly excised.

Stimulation Protocol
AF was induced by rapid pacing (5 V and a 125-190-ms interval) using bipolar electrodes implanted in the left atrium connected to a stimulation device (SEC-5104, Nihon Kohden, Japan). The duration of the stimulation pulse was set to 1 ms, the stimulation interval was started at about 190 ms and gradually shortened to about 120 ms (burst pacing), and the stimulation was stopped when the ECG waveform confirmed the induction of AF. The camera was manually triggered to start recording of optical mapping if AF was sustained after the stimulation was stopped ( Fig. 2 lower panel).

Cardiac Optical Mapping
Optical mapping was used to measure the cardiac excitation, which is the input signal for the conventional and proposed methods. To suppress motion artifacts caused by myocardial contraction, BDM (2,3-butanedione monoxime, 19 mM) was added to the perfusate. Next, the heart was then stained with 0.1-µM Di-4-ANEPPS, a voltage-sensitive dye. We used an LED ring light (primary wavelength: 520 nm, Xeon 3 Power Pure Green LED, OptoSupply, China) as the excitation source. Fluorescent signals passing through a long pass filter (passing > 600 nm, R-60, Nikon, Japan) were captured using a high-speed camera (1000 fps, 512-pixel × 512-pixel, FASTCAM Mini-AX50, Photron, Japan).

Analysis Method
We used the signals measured using the above experiment to compare the proposed and conventional methods. First, the membrane potential map was derived by inversion and normalization of camera signals at each pixel. To reduce noise, we applied a Gaussian filter with a size of 3 ms and 9 pixels in the temporal and spatial directions, respectively.

Fig. 2
Overview of the experimental system. Excitation of an excised porcine heart stained with a membrane potential dye is measured by capturing fluorescence using a high-speed camera. Bottom graph shows ECG signal during the induction of AF by applying burst pacing, and the subsequent manual camera trigger Subsequently, we identified the region of interest (i.e., 28.2 × 28.2 mm 2 ). Next, we applied both methods to the membrane potential map obtained as described above to detect FA. The details of the analysis for the proposed and conventional methods are the following. In the proposed method, the phase map is derived from the membrane potential map, which uses the Hilbert transform proposed by Bray et al. [7]. The derivation of phase values is as follows: Here, (t) is the phase value at time t and V′(t) is the membrane potential after baseline correction at each time point. Finally, H is the Hilbert transform. The resulting EWF map was also derived with θ act = − 2.73, σ = 0.54. After multiplying the EWF map by the divergence of the phase gradient, we binarized the images using different thresholds (0.7, 0.75, 0.8, 0.85, 0.9, and 0.95) and labeled the positive area. Finally, we calculated the center of gravity for each labeled area in the spatial and temporal directions to determine the occurrence position and time of FA.

Conventional Method
To evaluate the performance of the proposed method, conventional FA detection method based on local activation time (LAT) at each pixel was also applied to the same optical mapping recordings [16]. LAT of optical mapping at each pixel was determined by peak detection on the time derivative of the membrane potential signal. At each pixel and time, the time until the next activation was calculated to derive activation map (ACT). Then, CV vector was defined as the gradient of the ACT map and the divergence of CV was calculated.
Finally, we binarized the images using a threshold and labeled the positive area (i.e., as for the proposed method), calculated the center of the gravity of each labeled area in the spatial and temporal directions, and determined the position and time of FA occurrence.

Evaluation
A total of 126 FAs were visually observed during 9 episodes of 2-s optical mapping of the induced AFs. FA onset was intermittent and did not coincide with the stimulation .
frequency or the central frequency of the induced AFs. The accuracies of FA detection of the proposed and conventional methods of FA detection were compared. As the main evaluation metrics, precision, sensitivity, and F-measure are defined as follows: where TP is the true-positive rate, FP is the false-positive rate, and FN is the false-negative rate. Each detection result was determined to be a true positive if the spatial error between the detection position and the nearest true FA occurrence was < 4 pixels (4 mm) away and the temporal error was < 10 ms. When FAs were visually detected for accuracy evaluation, radial excitation propagation of at least about 3 frames had to be observed. Therefore, considering the propagation speed of WFs on optical mapping (about 1 pixel/ms), it was considered difficult to segregate FAs at distances less than 4 pixels, so the maximum allowable distance error was set to 3 pixels. In addition, since it is difficult to identify the exact timing of FAs by visual inspection, we set the maximum time error to 9 ms to account for the ambiguity of the correct timing. In case multiple detections for a single FA satisfied these conditions, any detection was counted as a true positive. Nine episodes (i.e., 9 × 2000 ms) of membrane potential signals measured via optical mapping were used as evaluation data. Due to the inability to detect the activation time before and after the limits of the measurement period, the activation map can fail to depict accurate excitation. To avoid this problem, we excluded the initial and final 200 ms of each episode from the evaluation data. Thus, we used a total of 14,400 ms of measurement signals for evaluation. Figure 3 shows an example of the comparison of the detection results from the conventional and proposed methods for the same FA. The images are the membrane potential maps measured by optical mapping (OPT), activation map (ACT) in the conventional method, phase map (Phase) and excitation wavefront map (EWF) in the proposed method, and the divergence map (Div) of the gradient of ACT and Phase in the conventional and proposed method, respectively. The white boxes enclose the position and timing of FA as , detected using each method. In the conventional method, multiple false-positive FAs were detected around the propagating excitation wave originated by both FA and planer wave from the left-top side. As supposed in Fig. 1, high divergence regions were observed on the WF lines and those resulted in false-positive errors. In contrast, the proposed method accurately detected a FA occurring in the center of the region and did not detect a false positive for the planar wave in the upper left region. In contrast to the conventional method, high divergence was not observed on the WF lines but at the center of FA in the proposed method, and the occurrence timing was precisely detected by multiplying EWF map to the divergence map of the phase gradient. Figure 4 shows an example of the detection results for each detection method against a planar wave. The left three images show membrane potential maps. The rightmost image displays the divergence of the CV vector field for the conventional method and the phase divergence (i.e., the multiplication result of the divergence of the phase gradient and the EWF map) for the proposed method. In the conventional method, we observed a high degree of divergence near the region where the WF passed through; accordingly, multiple FAs were detected incorrectly. In contrast, we also observed that the phase divergence calculated using the proposed method remained low even when the planar wave passed and false positives of FA events did not occur. Tables 1 and 2 show the quantitative evaluation results of the conventional and proposed methods, respectively. Figure 5 shows the precision and sensitivity of the conventional and proposed detection methods. For all binarization thresholds, the precision and sensitivity of the proposed method were higher than those of the conventional method. In the conventional method, the f-measure had a maximum value of 0.08 with a threshold value of 0.9 and precision and sensitivity were 0.05 and 0.26, respectively. On the other hand, the best precision and sensitivity values obtained via the proposed method were 0.46 and 0.91 at a threshold of 0.85 and 0.7, respectively. Moreover, the best f-measure value obtained using the proposed method was 0.58 at a threshold of 0.8.

Effects of the Proposed Method on FA Detection
As shown in Fig. 5, the precision and sensitivity of the proposed method were greater than the conventional method. These results suggest that the proposed method has the potential to suppress false positives of FA readings that are observed in the conventional method and thereby to improve detection accuracy. In the following section, we discuss the reasons for the improved precision and sensitivity of the proposed method, respectively.
In the conventional method compared in this study, a CV map of the field of view was created from the LAT derived from the optical signal. However, to obtain a CV map using the conventional method without the discontinuity problem shown in Fig. 1, the analysis period and ROI should be selected appropriately so that the ROI does not contain a WF in the beginning and end frames of the period and that the WF passes all the way through the ROI in the analysis period, as shown in Fig. 6. On the other hand, a method has recently been proposed that does not rely on LAT and calculates CV using the correlation of signals at multiple measurement locations. However, because complex AFs are inherently non-stationary and not repeatable, the same analysis period and ROI limitations as described above are required to analyze CVs, which are updated with each WF passage.
These limitations in period and ROI significantly facilitate the automatic detection of FAs, but it is difficult to select periods and ROIs that satisfy the above conditions during a complex AF, such as the one evaluated in this study. Thus, although the conventional method can be said to quantify the indicators representing FAs, it cannot be said to fully automate the detection of FAs in complex AFs. As a result of this experiment, which was conducted without such period and ROI limitations, we believe that we have clarified the inherent performance limitations of the conventional method in the automatic FA detection performance during complex AFs and have demonstrated the superiority of the proposed method.
We attribute the significant improvement in the precision of the proposed method to the use of the phase value as a circular variable to detect FA. This helps to avoid inherent problems associated with the discontinuity near the WF of the activation map, which is a fundamental limitation of the conventional method. Figure 4 shows that in the conventional method, the divergence value was high along the WF of the planar wave and resulted in multiple false-positive FAs. This is because, as described in Sect. 2, the direction of excitation propagation could not be correctly evaluated from the spatial gradient of the activation map due to the discontinuous boundary near the WF. In contrast, the proposed method correctly evaluated the direction of excitation propagation from the spatial gradient of the phase value (i.e., the circular variable), and the phase divergence value that would have otherwise existed was suppressed. Therefore, the proposed method, which uses the phase value as a circular variable, suppressed false positives of FA observed in the conventional method and remarkably improves the precision of FA detection. Fig. 4 An example of the detection of a planar wave using each method. The left three images are membrane potential maps. The rightmost image shows the divergence of the CV vector field for the conventional method and the phase divergence (i.e., the multiplication result of the divergence of the phase gradient and EWF map) for the proposed method  In addition, the proposed method also improved the sensitivity of detection relative to the conventional method. Missed detection of FAs by the conventional method occurs when the method fails to detect passing excitation at the FA occurrence position on the activation map (supplementary figure). In the derivation of the activation map in this study, the activation time was determined to be the peak of the membrane potential change rate dV/dt, as described in Sect. 3.4. Moreover, it is known that the resting potential of myocardial tissues can increase during periods of complex excitation, such as AF, and the amplitude of the membrane potential at this time will be lower [17]. This can lead to oversight of the excitation passage by the conventional method. However, in the proposed method, the phase value is calculated from the membrane potential using Eq. 3, and the independence of the excitation passage detection from the rate dV/dt contributes to the improved sensitivity.
As our discussion shows, the improved sensitivity observed in this study is not an essential effect of the proposed method using the phase value as a circular variable. Rather the improvement relates to the experimental method of this study, which uses the membrane potential, a clinically difficult measurement. However, decreasing amplitude of membrane potential signals during AF can lead a detection protocol using the conventional method to miss excitation events and are therefore a common feature observed in clinically measurable unipolar and bipolar electrode signals [18]. When applying the proposed method in clinical settings, it is necessary to derive the phase map using a robust method for detecting excitation passage, even for low-amplitude electrode signals during AF.

Analysis of Error Factors
The detection accuracy of FA improved when using the proposed method. However, some false positives and false negatives remained unresolved. With a binarization threshold of 0.8, the proposed method yielded 72 cases of false detection and 20 cases of missed detection. One case of missed detection was due to an error in the matching of detection results with the true results, as described in Sect. 3.5. However, the remaining 91 cases were all errors related to FA events at the WF of the planar wave, as described below. Figure 7 shows an example of a local FA located near the WF of the planar wave observed in the evaluation experiment. In this example, an increase in membrane potential was observed near the WF of a planar wave propagating from the upper left direction around 0 ms. Moreover, there was also a slight radial excitation observed between 0 and 3 ms. At 4 ms, the excitation merged with the planar wave. All 72 cases of false positives observed in this study were caused by the local FA near the WF of this planar wave and were missed during the visual labeling of the true results. In contrast, in 19 cases of missed detections, visual labeling identified the local FA near the WF. However, conduction was observed numerically between the WF and the occurrence position of FA and led to missed detection errors.
In summary, the main cause of error with respect to missed detection was the ambiguity inherent in manually labeling occurrences of localized and momentary FA events near the WF of a planar wave. Moreover, the results of the proposed quantitative method are more reliable than those that require manual labeling. The presence of such FA events has not been previously recognized, and the importance of the excitation revealed by the proposed method enabled the quantitative detection of FA events during complex excitation periods that were precisely measured via optical mapping. However, the mechanism behind the observed local FA events remains unclear. Possible hypotheses related to the underlying mechanism include the possibility of wraparound of three-dimensional scroll waves due to fast passage deep within tissue or heterogeneous depolarization thresholds in response to electrotonic effects associated with changes in the intracellular state during fibrillation. Further research is necessary to clarify the underlying mechanisms. In addition, it is necessary to determine whether such FA events can be effective targets for ablation therapy. If so, a high-resolution electrode map would be necessary to distinguish between the planar wave and nearby FA events.

Limitations
In this study, we used optical mapping, which enables measurement of the membrane potential with high temporal and spatial resolution, for the evaluation of the proposed method. However, optical mapping cannot be used in clinical settings because it uses toxic dyes. When using the proposed method in clinical settings, it is therefore necessary to derive the phase map from electrode signals that are easier and safer to obtain in clinical settings. However, electrode measurements typically have lower spatial resolution than optical mapping. When applying the proposed method to phase maps derived from electrode signals, there is a possibility of missing FA events that occur between electrodes. Moreover, there is also a possibility that detection errors by the Fig. 6 Required conditions of the analysis period and ROI for conventional methods. a Scheme of the required ROI. Top: a case of one-way planer wave. Bottom: a case of focal activation. In both cases, it is required that the ROI does not contain a WF in the beginning and end frames of the period and that the WF passes all the way through the ROI in the analysis period. b Example of the required ROI during complex fibrillations. While the ROI of rectangles with white solid line satisfies the conditions, the ROI including rectangles with dashed lines does not satisfy Fig. 7 An example of a local FA near the WF of a planar wave. An increase in membrane potential was observed near the WF of a planar wave propagating from the upper left direction at ~ 0 ms. A slight radial excitation was observed at ~ 3 ms, and it then merged with the planar wave at 4 ms proposed method can be attributed to hum noise derived from AC power sources. In contrast, in recent years, mapping techniques have been developed to visualize excitations with high accuracy using electrode signals [9,19,20]. By combining such mapping techniques with the proposed method, it may be possible to accurately detect the position and time of FA events even in clinical settings.

Conclusion
In this study, we proposed a quantitative detection method for FA using the divergence of the spatial gradient of the phase map, which was a circular variable. We evaluated the detection accuracy of the proposed method using an excised porcine heart. Our evaluation results indicated that the precision and sensitivity of the proposed method were considerably higher than the conventional method. The proposed method showed potential to suppress false positives of FA events that are observed in the conventional method.