Anticipation of ventricular tachyarrhythmias by a novel mathematical method: Further insights towards an early warning system in implantable cardioverter defibrillators

Implantable cardioverter defibrillators (ICD) are the most effective therapy to terminate malignant ventricular arrhythmias (VA) and therefore to prevent sudden cardiac death. Until today, there is no way to predict the onset of such VA. Our aim was to develop a mathematical model that could predict VA in a timely fashion. We analyzed the time series of R-R intervals from 3 groups. Two groups from the Spontaneous Ventricular Tachyarrhythmia Database (v 1.0) were analyzed from a set of 81 pairs of R-R interval time series records from patients, each pair containing one record before the VT episode (Dataset 1A) and one control record which was obtained during the follow up visit (Dataset 1B). A third data set was composed of the R-R interval time series of 54 subjects without a significant arrhythmia heart disease (Dataset 2). We developed a new method to transform a time series into a network for its analysis, the ε−regular graphs. This novel approach transforms a time series into a network which is sensitive to the quantitative properties of the time series, it has a single parameter (ε) to be adjusted, and it can trace long-range correlations. This procedure allows to use graph theory to extract the dynamics of any time series. The average of the difference between the VT and the control record graph degree of each patient, at each time window, reached a global minimum value of −2.12 followed by a drastic increase of the average graph until reaching a local maximum of 5.59. The global minimum and the following local maxima occur at the windows 276 and 393, respectively. This change in the connectivity of the graphs distinguishes two distinct dynamics occurring during the VA, while the states in between the 276 and 393, determine a transitional state. We propose this change in the dynamic of the R-R intervals as a measurable and detectable “early warning” of the VT event, occurring an average of 514.625 seconds (8:30 minutes) before the onset of the VT episode. It is feasible to detect retrospectively early warnings of the VA episode using their corresponding ε−regular graphs, with an average of 8:30 minutes before the ICD terminates the VA event.


Introduction
Implantable cardioverter defibrillators (ICD) are the cornerstone of sudden cardiac death prevention through termination of ventricular tachycardia/ventricular fibrillation. Although ICD shocks usually occur when the subject is unconscious, it could be very useful to patients and close relatives to have the possibility to know in advance, either seconds or minutes, when those malignant arrhythmias could occur in order to take appropriate preventive measures. We hypothesized that a novel mathematical analysis, ε−regular graphs, could perform such task.
Network theory possesses the capacity to abstractly represent interactions of any kind of entities. Currently, complex networks have arisen as a common way to tackle intricate dynamics [1]. A broad range of applications in different biological and medical areas abound. In the area of biology, they have been used to analyze a population's structure [2,3], and pandemics [4,5]. The use of protein-protein interaction networks coupled with information theory have led to discover potential therapeutic biomarkers on cancer research [6]. Integrative approaches for anticipating critical transitions have been proposed [7] in several phenomena, although the area of cardiology has not yet been explored. The terms "early warnings" and "tipping points" are still not part of the cardiologist community. Several methods have been developed to transform time series into networks for its analysis. Such methods include the visibility graphs method [8], and a plethora of its modifications [9,10], which consider the topological properties of the time series, the recurrence analysis of time series [11,12], and the analysis based on the phase space [13]. In this work, we usher in a new method to transform a time series into a network for its analysis, the ε−regular graphs. This novel approach transforms a time series into a network which is sensitive to the quantitative properties of the time series, it has a single parameter (ε) to be adjusted, and it can capture long-range correlations. This procedure permits using graph theory to extract the dynamics of any time series. As a direct application of ε −regular graphs, data from patients diagnosed with imminent ventricular tachyarrhythmias (VT) was analyzed. The heart activity is driven by the action of the opposing forces of the sympathetic and the parasympathetic nervous systems [14,15]. The failure in heart function is the result of malfunctions in the myocardium, heart valves, pericardium, or the endocardium [16].

Limitations
The method proposed in this work requires further testing with patients whose clinical history is well documented and controlled, coupled with respiratory data. A significant clinical limitation of this work is the fact that this approach is restricted to patients with normal sinus rhythm and is unlikely to work in patients with atrial fibrillation or those with a pacemaker. In the former group because of the large variability of R-R intervals, and in the later because of fixed pacing rhythms.

A mathematical method to transform time series to networks
The method consists in assigning to each point in a time series a vertex in the network. Then, for a fixed value of the parameter ε = ε 0 , any two points of the time series p 1 ,p 2 will be joined in the network, if and only if |p 1 −p 2 |�ε 0 ; this means that two point of the time series will be adjacent in the ε−regular graph if the values of the points have a maximum difference of ε 0 . For illustrative purposes, we show a diagram of the algorithm in Fig 1. In Fig 1A, we show a time series of ten points; the values of the points are p 1 ,p 7 , and p 10 = 0.2; p 2 and p 4 are equal to 0.29; p 5 = p 8 = 0.38; p 3 = p 6 = p 9 = 0.7. The ε−regular graph is constructed with a parameter Tachyarrhythmia Database Version 1.0 from Medtronic Inc. (available at doi.org/10.13026/ C25K5D). A set of 81 pairs of R-R interval time series records from different patients was obtained (Group 1A), each pair contains one record before the VT episode and one control record (CR) which was obtained during the follow up visit (Group 1B). A third data set composed of the R-R interval time series of 54 subjects without a significant arrhythmia heart disease was obtained (available at doi.org/10.13026/C2NK5R). The authors did not have special access privileges.

Competing interests:
The authors have declared that no competing interests exist.
value of ε 0 = 0.1. In Fig 1B, intervals of width ε 0 = 0.1 are drawn around each point of the time series. For a given point p of the time series, all the point lying inside the interval of width ε 0 will be adjacent to p in the corresponding ε−regular graph.
Other algorithms to convert time series into graphs have been developed but with qualitative instead of quantitative rules for determining the adjacencies in the corresponding graphs (visibility plots). The algorithm for the visibility graphs determines the adjacencies by analyzing the lines joining the points of the time series [8] as observed in Fig 1C. The visibility graph algorithm confers to its graph properties that strongly differ to our ε−regular graph derived from the same time series because the former does not capture the quantitative properties of a width ε 0 will be adjacent to p in the corresponding ε−regular graph. Note that points with a higher value of 0.7, belong to a different component than the rest of time series reflecting its outlier nature. Also note that the periodic values, p 3 ,p 6 and p 9 form a subgraph. Points with the same value that are not periodic will form a complete subgraph, such as the points p 1 ,p 7 and p 10 . In (C), the time series is transformed to a graph using the visibility-graph algorithm.
https://doi.org/10.1371/journal.pone.0235101.g001 time series as the ε−regular graph do. When considering the outliers of a time series, in a visibility graph, the extremely high or low values of a time series would be "visible" from almost all the rest of the points of the time series, and thus, according to the visibility graph algorithm would be highly connected and act as a hub. This would remain true even if the outlier value were not so extremely contrasting. In contrast, an outlier in a ε−regular graph would have different behaviors according to its value. If the outlier value is significantly different, it will be assigned to a different component in the ε−regular graph: If the suspected outlier's point value is not so different to the rest of the time series, it will remain in the same component. This behavior is reflected in the time series of Fig 1, where points with a higher value of 0.7, belong to a different component than the rest of time series reflecting its outlier nature. If instead of the value 0.7, the values were set to 0.48, they will still be higher than the rest of the time series but will remain in the same graph component of the rest of the time series points since the ε 0 value is set at ε 0 = 0.1, and the points with the second-highest value are equal to 0.38.
From the definition of adjacencies in an ε−regular graph, it is directly derived that on a given set of points, in which the points of the set have a value difference up to ε 0 among them, they will be joined and thus they will form a complete subgraph. This result is useful when considering time series with regular or periodic values. The periodic values will form a complete subgraph, as the points p 3 ,p 6 , and p 9 in Fig 1B. Also, points with the same value that are not periodic will form a complete subgraph, such as the points p 1 ,p 7 , and p 10 (Fig 1) that form the complete graph K 3 . A similar property is not inherited in visibility graphs. In visibility graphs, periodic or regular points from the time series might in some cases not be adjacent, as some points might block the visibility condition for them to be adjacent. What can be rescued from visibility graphs is the short-term correlations of the time series.
As proof of concept, two time series were simulated from theoretical frameworks and their corresponding ε−regular graphs were derived. Time series of 10 4 points were simulated, the first time series was obtained from a standard normal probability distribution, and the second from a standard Brownian motion valued at integer times. The value of the parameter ε 0 was set as the standard deviation, and half of the standard deviation of both the normal distribution and the Brownian motion, which correspond to 1 and 1/2 respectively, in both cases. The degree centrality measure was calculated for the resulting graphs (Fig 2). The values of the centrality measure were standardized to the unit interval. The degree distributions for the graphs constructed form the Brownian motion approximates a Gaussian curve despite its multiple modes; and the distribution from points sampled from a normal distribution approximates a lognormal distribution. The distributions of the degree measurements maintain, up to some extent, the statistical properties from the time series they were derived from. The difference between any two values of the Brownian motion follows a normal distribution, and this property is shown by the ε−regular graphs in the degree centrality measure when setting the ε−parameter equal to values related to the standard deviation.
The computer program is found in Supporting Information I.

Application to imminent VT
As a direct application of the ε−regular graphs algorithm, a set of time series related to VT was analyzed. The data was downloaded from the Spontaneous Ventricular Tachyarrhythmia Database Version 1.0 from Medtronic Inc. (available at http://physionet.org/physiobank/ database/mvtdb/) [17]. A set of 81 pairs of R-R interval time series records from different patients was obtained (Group 1A), each pair contains one record before the VT episode and one control record (CR) which was obtained during the follow up visit (Group 1B). A third data set composed of the R-R interval time series of 54 subjects without a significant arrhythmia heart disease was obtained (available at: https://physionet.org/physiobank/ database/nsrdb/) [16]. This third dataset of time series will be hereafter considered as healthy subjects (HS) and denoted as Group 2.
The time series for groups 1A and 1B were cropped to the same length starting from the end to make them directly comparable (985 points), the start of the VT episode is at the last point of the time series. The Group 2 (HS) series were also cropped by subsampling a random set of sequential points of each subject that match the length of the VT and CR time series. The VT, CR, and HS time series of each patient were subdivided in time series of 60 points with a sliding window method with an offset of one point. This would allow analyzing the change of the series in time to detect an early warning of the VT episode. The time series on each window were transformed into graphs with the ε−regular graph algorithm setting the parameter value at ε 0 = 0.04.
The average degree of the graphs from each window of each subject was calculated and averaged among the subjects, for the three different time series. This process results in a time series of the average degree of the graphs representing the three different states of the subjects (Fig 3). The datasets VT and CR arise from the same subjects, so a direct comparison of subjects prior to a VT episode and in a normal stage is possible. The average of the difference between the VT and the CR graph degree of each patient, at each time window (Fig 3), reaches a global minimum value of −2.12, followed by a drastic increase of the average graph until reaching a local maximum of 5.59. The global minimum and the following local maxima occur at the windows 276 and 393, respectively. This change in the connectivity of the graphs distinguishes two distinct dynamics occurring during the ventricular tachyarrhythmia, while the states in between the 276 and 393, determine a transitional state. We propose this change in Fig 2. Degree centrality of the applied ε−graph algorithm. Two time series were simulated from theoretical frameworks and their corresponding ε-regular graphs were derived. Time series of 10,000 points were simulated, the first time series was obtained from a standard normal probability distribution (solid and dashed yellow curves), and the second from a standard Brownian motion (solid and dashed blue curves) valued at integer times. The value of the parameter ε 0 was set as the standard deviation, and half of the standard deviation of both the normal distribution and the Brownian motion, which correspond to 1 and 1/2, respectively, in both cases. The degree distributions for the graphs constructed form the Brownian motion approximates a Gaussian curve despite its multiple modes; and the distribution from points sampled from a normal distribution approximates a lognormal distribution. The optimization of the parameter ε is based on the statistical parameters of the R-R time series. The average, minimum and maximum distance between two points of the R-R time series are: 0.043, 1.11×10 −16 , 0.12, respectively. Thus, the ε value of 0.04 approximates the mean of the differences. In Supporting Information II, we show the degree time series when the values 0.01, 0.02, 0.03, 0.04, and 0.06 are assigned to the parameter ε (S1 Fig). Note that the overall pattern is preserved. In particular, the abrupt change of the dynamics is captured by the different values of the parameter ε.

Comparison with detrending fluctuation analysis
A widely used procedure used in the analysis of data originated from diverse heart records is the application of the Detrended fluctuation analysis (DFA) method. DFA is a mathematical linear method to analyze time series by removing the linear trend of time series divided into smaller windows. This method is of special use to address nonstationary time series [18]. Results from the DFA method are commonly graphed in a log-log scale and the scaling exponent of the time series is estimated from a least-square fit of a linear model. The scaling exponent measures the correlation in the noise and approximates the Hurst exponent of the time series.
The DFA method was applied to the mean time series of the 3 groups (Fig 3A), and the difference between the mean VT and CR subjects (Fig 3B). From the DFA of the different states of the subjects it can be observed that the time series from the HS subjects possess the same scaling properties at short-and long-time lengths, which is deduced from the fact that the DFA approximates a linear model. On the other hand, the DFA analysis from the VT and CR subjects show that their time series possess two different scales of autocorrelation, which is related to the two different linear models fitting the DFA of VT and CR. The VT and CR time series behave similarly in the sense that both exhibit different scaling properties for short correlations (windows with 15 or less points) and a different scaling factor for long correlations (15 or more points). Since the slope of the liner models fitted to VT, 0.24 and 0.95 for short and long correlations, respectively, and CR, 0.25 for short and 1.01 for long correlations ( Fig  3C), are practically the same, then, their corresponding Hurst exponents will be the same, which results in that the VT and CR behave similarly regardless if short or long correlations are assessed. This is validated by the fact that the slopes of the two linear models fitted to the DFA of VT-CR, 0.31 and 1.08 for short and long correlations, approximates the ones obtained when analyzing VT and CR separately (Fig 3D). The DFA method is capable to discern the different scaling properties occurring on the Groups 1A (VT) and 1B (CR) patients as compared to the Group 2 (HS) subjects. However, the DFA results are not varying in time, and hence this method is not capable of discerning an early warning for VT.

Discussion
In this work we propose a novel parametric method to analyze time series by transforming them into networks. By using this method, it is possible to apply the graph and network theory in the analysis of time series. Herein, a direct application of the ε−regular graph method is herein shown by using time series data derived from patients with ventricular heart tachyarrhythmia disease. The application of the ε−regular graph method, using a sliding window framework, detected a potential early warning of the disease that it is not detectable using the current linear methods available for the analysis of time series. The ε−regular graphs differ from the visibility graph method as the former is a parametric quantitative method and the latter is a qualitative approach. The adjustable parameter ε in ε−regular graphs, determines the sensitivity of the transformation of time series to networks. By varying the parameter, it is possible to obtain a range of graphs going from graphs in which a vertex is only connected to other ones having the same value, up to completely connected graphs. An inverse transformation, form a network to a time series, would be possible if there exists a compendium of graphs derived from the same time series using different ε values. Then, if needed, the original time series can be inferred using the different adjacencies from the ε−regular graphs. An inverse transformation that faithfully recovers the time series is not possible for visibility graphs. Since visibility graphs is a qualitative methodology, the values of a time series derived from these graphs would vary in an interval, whose length would be different for each point in the time series. The framework of complex networks for analyzing heart rate variability data towards the detection of early warnings and the design of clinical tools for disease management has been considered before as other nonlinear methods [19]. Visibility graphs have been applied to the analysis of congestive heart failure [20]. Inhere, a statistical analysis of the scale-freeness of the obtained network is used for the detection of early stages of the disease. In a broader analysis, several summary statistics of a horizonal visibility network have been proposed as useful for the analysis of heart rate variability [21]. In general, the use of summary statistics for the detection of early warnings in a transition of dynamical state may be difficult since such statistics may rely on inadequate data or other factors [22]. Other studies have shown that the incorporation of respiration signals to the electrocardiogram data increase the detection of a VT episode [23]. Hitherto, the effect of the vagus nerve in the heart activity has been recently investigated [25]. Different techniques based on other methodologies and data have shown different times before the VT episode occurs [23,24]. Any predictor, regardless of the methodology must clearly distinguish a VT episode from the usual cardiac arrythmias of each patient to avoid false positive detection. So far, the low heart rate variability has been considered as the single predictor of heart failure, although the forces for the acceleration and deceleration in heart activity have been shown to be uncoupled [25]. The device used by the patient has high impact on any method for the detection of early warnings of a cardiac malfunction, as it has been shown that ICDs can detect QT variability in near-field or far-field right ventricular intracardiac electrogram [26]. ICD are excellent machines devoted to terminating VT and they have proved its efficacy to prevent sudden cardiac death in different clinical settings. The performance of the algorithms has been tested first to detect the VT and to provide appropriate shocks. Then, algorithms were improved to avoid unnecessary ("inappropriate") discharges to the patient. In recent years there has been a small but strong movement in the medical community towards the possibility of alerting the patient when an ICD shock is going to occur. This possibility is not minor. From a clinical point of view, such alert could permit the patient or his close relatives to take appropriate measures before the shock takes place. A new window of opportunity (clinical interventions) could be generated if a software could be able to detect with some seconds or, even better, minutes, the possibility of an imminent ICD shock. Until today, there is no such possibility. The present retrospective study sheds light of a possible mathematical analysis that could detect "early warnings" of an appropriate ICD shock for VT with an average of 8:30 minutes. The process of optimization of the ε parameter value requires a more extensive clinical experimentation. It stands to reason that the parameter value is specific for each individual patient and it ought to be tuned from the complete set of clinical parameters of the patient to avoid false-positives and false-negatives. In a more general case, it is also probable that the parameter value of ε is not fixed throughout the day and is dependent of the circadian rhythm of the patient.
Obviously, this mathematical application should be tested prospectively, but this can only be done if implemented into the software of the ICD. Collaboration with ICD industry is vital to achieve such goal.

Conclusions
Early warnings of the VA episode could be detected using their corresponding ε−regular graphs, even 8:30 minutes before the ICD comes into action. A prospective study is warranted to further corroborate this finding.