Bimodality and Scaling in Recurrence Networks from ECG data

TL;DR

This study uses recurrence network analysis on ECG data to reveal bimodality in degree distribution and scaling behaviors, providing insights into cardiac dynamics and potential disease markers.

Contribution

It introduces the first detailed analysis showing bimodality and scaling in ECG recurrence networks, linking these features to complex cardiac dynamics and disease variations.

Findings

Recurrence networks from ECG exhibit bimodal degree distributions.

Two distinct scaling regions are identified in link density.

Disease-specific variations in bimodality and scaling are quantifiable.

Abstract

Human heart is a complex system that can be studied using its electrical activity recorded as Electrocardiogram (ECG). Any variations or anomalies in the ECG can indicate abnormalities in the cardiac dynamics. In this work, we present a detailed analysis of ECG data using the framework of recurrence network (RN). We show how the measures of the recurrence networks constructed from ECG data sets, can quantify the complexity and variability underlying the data. Our study shows for the first time that the RN from ECG show the unique feature of bimodality in their degree distribution. We relate this to the complex dynamics underlying the cardiac system, with structures at two spatial scales. We also show that that there is relevant information to be extracted from the scaling of measures with recurrence threshold. Thus we observe two scaling regions in the link density for ECG data which is compared with scaling in RNs from standard chaotic and hyperchaotic systems and noise. While both bimodality and scaling are common features of RNs from all types of ECG data, we find disease specific variations in them can be quantified.

Tables

Table 1. Table 1 : Network measures, CC, APL and LD, for ECG data sets, standard nonlinear systems and white noise.

Data	CC	APL	LD
Healthy	0.825 $\pm$ 0.034	4.439 $\pm$ 0.689	0.338 $\pm$ 0.103
BB	0.741 $\pm$ 0.071	5.364 $\pm$ 2.507	0.199 $\pm$ 0.157
CM	0.821 $\pm$ 0.093	4.281 $\pm$ 1.688	0.390 $\pm$ 0.222
DR	0.837 $\pm$ 0.044	4.167 $\pm$ 1.134	0.411 $\pm$ 0.119
MI	0.811 $\pm$ 0.066	4.255 $\pm$ 0.973	0.319 $\pm$ 0.164
Chen	0.573 $\pm$ 0.005	5.611 $\pm$ 0.266	0.018 $\pm$ 0.002
R$\ddot{o}$ssler	0.624 $\pm$ 0.003	8.985 $\pm$ 0.071	0.018 $\pm$ 0.000
White Noise	0.477 $\pm$ 0.008	3.213 $\pm$ 0.126	0.034 $\pm$ 0.005

Table 2. Table 2 : Scaling indices for link density calculated from log -log plots for ECG data of all types, standard nonlinear systems and white noise.

Data	Scaling index ${\gamma }_1$	Scaling index ${\gamma }_2$
Healthy	0.782 $\pm$ 0.288	0.387 $\pm$ 0.296
BB	1.248 $\pm$ 0.509	0.830 $\pm$ 0.364
CM	0.891 $\pm$ 0.621	0.407 $\pm$ 0.401
DR	0.609 $\pm$ 0.271	0.426 $\pm$ 0.368
MI	0.941 $\pm$ 0.449	0.372 $\pm$ 0.276
Chen	2.462 $\pm$ 0.059	2.136 $\pm$ 0.199
R$\ddot{o}$ssler	1.859 $\pm$ 0.010	1.263 $\pm$ 0.017
White Noise	1.451 $\pm$ 0.157	0.000 $\pm$ 0.000

Equations

$$LD=(\frac{1}{N^2})\sum _{i,j=1}^N{A}_{ij}$$

$$C_i=\frac{{\sum }_{j,q}{A}_{ij}{A}_{jq}{A}_{qi}}{k_i(k_i-1)};CC=\frac{{\sum }_{i=1}^N{C}_i}{N}$$

$$APL=\frac{1}{N(N-1)}\sum _{i\ne j}{l}_{ij}$$

Full text

1 Introduction

The complex systems occuring in nature are mostly studied using observational or measured data over time and the methods of nonlinear time series analysis are very useful in analysing the data to understand their complexity[1]. However, for short data sets, recent technique of recurrence networks is found very effective[2, 3]. In this study we use the approach of recurrence networks to study the complexity of heart dynamics using the clinical Electrocardiogram (ECG) data of one minute duration. Such a study can bring out variations or anomalies in the ECG that carry information about the abnormalities in the underlying cardiac dynamics.

Recently attempts have been made to understand ECG signals from a dynamical systems perspective establishing low dimensional chaos in the cardiovascular system [4]. Also, significant advances have been made in the direction of cardiac modelling[5, 6] which can generate artificial signals mimicking features of a real ECG. Most of the studies reported in this direction are on heart rate variability (HRV) data to classify cardiac abnormalities[7, 8] and the analysis of the full ECG waveforms remain relatively under-explored[9]. Using multifractal measures it is reported that, the variability in the complexity of the cardiac dynamics is less in the case of healthy subjects as compared to patients[10, 11]. But,the conventional techniques of nonlinear time series analysis for multifractal measures require long data sets and hence not very useful for short duration clinical ECG data.

In this study we report the interesting features that we find in the characteristics of Recurrence Networks (RNs) from short duration clinical ECG data. RN, by construction, reflects the recurrence pattern of points of the attractor or trajectory in the phase space, reconstructed from the data. Hence the measures characterising the network such as degree distribution, clustering coefficient, average path length, link density etc. can effectively quantify the complexity and variability in the ECG signals.

Our study includes mainly ECG data from 96 cases selected from the PTB diagnostic database[12] of Physionet[13], within the age group of 30-75 years. Among these 33 are from healthy volunteers while the rest are classified into one of the following cases: Bundle Branch Block (BB), Cardiomyopathy (CM), Dysrhythmia (DR) and Myocardial Infarction (MI). We construct RN from each ECG time series and study how the network measures differ for typical cases. We also study how the measures vary when the threshold used in constructing the recurrence network is varied. This can provide a finer level of characterisation of the reconstructed phase space trajectory and hence the nature of the underlying dynamics. We compare our main results from ECG data sets with similar measures derived from data of standard chaotic and hyperchaotic systems and white noise.

Our results for the first time bring out a novel feature of the RNs from ECG data, viz. the bimodal nature of its degree distribution. We also observe a typical scaling behaviour for link density as a function of the recurrence threshold, indicating a specific structure for the phase space attractor. While both these features are observed in general for RNs from all types of ECG data, we show their quantifiers can be useful as distinguishing measures to detect disease cases from healthy. This is because they can capture the information about any abnormalities or variations in the complexity of the underlying dynamics among the datasets.

2 Construction of Recurrence Networks from embedded phase space attractors of ECG data

The first step in our analysis is to reconstruct the phase space attractor from the ECG data using the method of time delay embedding[14]. As is well known, the ECG recordings are prone to noise and artefacts due to body movement, power-line interference, baseline wander and muscle contraction etc[15]. So we do pre-processing of data to remove artefacts and trends by passing them through a band pass filter (0.5-50 Hz)[16]. The filtered time series is then normalized to the interval [0,1] [10]. The data thus obtained has 60000 points, which is down-sampled to 6000 points ($x_{t}$) by binning. We make sure that the whole process retains all the essential features of the original data sets.

Following conventional approaches, we take the time delay for embedding as the value at which the autocorrelation of data points, falls to $\frac{1}{e}$ of its initial value. The embedding dimension $m$ is chosen using the method of False Nearest Neighbours (FNN)[17]. We find for most of the data sets used, the FNN gives $m=4$ or less and hence for uniformity we embed all data in a space of dimension $m=4$. With the time delay $\tau$ and embedding dimension $m$ thus fixed, we get the vectors in the reconstructed phase space for each data[18]. A 2-dimensional projection of the reconstructed trajectory or attractor in phase space for typical cases of ECG data are shown in Fig. 1. For healthy ECG, we observe a structure with large and small loops. In the case of BB and DR, the attractors are much wider while in the case of CM, the small scale structure is much reduced compared to healthy. The attractor has a richer structure in case of MI. We discuss how these features relate to the network measures in the following sections.

Among the various methods reported to generate complex networks from a given time series[3], we find the recurrence networks[2] is the most suitable for the short-duration 1-minute ECG data. The RN is constructed to capture the pattern of recurrences of points or dynamical states in the reconstructed phase space trajectory. This is done by defining each point as node of the network and points that recur within a chosen threshold $\varepsilon$ are connected by links. The occurrence of links is represented in the form of a binary matrix called Adjacency matrix defined as $\displaystyle A_{ij}=R_{ij}-\delta_{ij}$, where $\displaystyle R_{ij}=\Theta(\varepsilon-\left\|\vec{v_{i}}-\vec{v_{j}}\right\|)$ and $\Theta$ is the Heaviside step function[2]. All the required measures representing the complexity of the network are then computed from the adjacency matrix $A$.

The choice of $\varepsilon$ has been discussed extensively in the literature[3]. The rationale behind all of them is that the resulting network should reflect the important recurrences, but should not result in an over-connected network that will mask the relevant features. We observe that for ECG time series considered, most of the networks become just connected at $\varepsilon$ = 0.1, making it a suitable choice across all data sets. The RNs thus constructed for typical ECG data sets are shown in Fig. 2.

We mainly study how to characterise the RNs using degree distribution, link density, clustering coefficient and average path length[19]. The degree of a node is defined to be the total number of its connections or links. For a given node $i$, we calculate degree $k_{i}$ as $k_{i}=\sum_{j=1}^{N}A_{ij}$, where $N$ is the number of nodes in the network and the degree distribution $P_{k}$, is the probability of nodes with a particular degree $k$.

The link density (LD) is computed as the ratio of actual number of links in the network versus total number of possible links using,

[TABLE]

The clustering coefficient (CC) for the network is defined to be the average of the local clustering coefficients. The local clustering coefficient for a node $i$ is defined to be the ratio of the number of triads it is part of, and the number of such possible triads including the node $i$.

[TABLE]

It is clear that every pair of nodes in the network can have one or more possible paths and the shortest path between a pair of nodes (i,j) corresponds to the minimum number of links connecting them. The average path length (APL) for the network is defined as the average over the shortest distances for every pair of nodes in the networks.

[TABLE]

where $l_{ij}$ is the shortest distance (minimum number of links) between nodes $i$ and $j$.

In the following sections, we discuss how these network measures capture the underlying dynamics of the cardiac system.

3 Measures of Recurrence networks from ECG data

3.1 Degree distribution

The simplest measure for the structure of connections in a network is the degree distribution and complex networks are often classified using this single measure. Thus scale-free networks from real systems exhibit power law degree distribution and random networks mostly have Poisson distribution[19].

For RNs from ECG data, we observe a novel feature, bimodality in the distribution for degrees and to the best of our knowledge this has not been reported so far for recurrence networks derived from time series. The results presented are from channel-2 data from each ECG and we find that the other channels also lead to qualitatively similar results. The presence of two peaks with almost negligible points in between seems to be a common feature of RNs from all types of ECG data. But the distribution appears to be distorted for data from unhealthy cases compared to healthy. For comparison we also show the degree distributions from RNs of standard nonlinear dynamical systems, Rössler system (chaotic) and the Chen system (hyperchaotic)[20], and that from white noise, which are all unimodal in nature.

For RNs from healthy data, the bimodal nature can be correlated to the phase space structure in Fig. 1, where we note the presence of a small dense loop like region and a large ring. The presence of the dense region is what gives rise to the second peak at higher degrees. In the case of CM, the attractor has a highly dense small region with a large ring. For BB, the distribution is unimodal at this threshold. This must be due to the absence of small structure and a wider ring structure in phase space. The degree distribution for DR is bimodal. For MI, the first peak is broad and the distribution is bimodal.

In comparison, for white noise, the distribution is a single wide peak, indicating uniform filling of the attractor, without any small scale structures. In the case of a standard chaotic system like Rössler system, we observe the average degree to be much lower than that of the noise and the distribution is unimodal. This behaviour can be understood in terms of the distribution of phase space points confined to a smaller region in phase space. For Chen system, in hyperchaotic region, we observe a fat tailed distribution, with a wider range of degrees but unimodal character which indicates absence of a well defined two-scaled structure in the attractor, in contrast to the ECG. Thus we see that the finer structures and complexity of the attractor is reflected in the degree distribution and in that respect, the RNs from ECG form a different class of networks.

We compute values of CC, APL and LD from RNs of all the data sets and present the results in Table 1. We observe that in general, the RNs from ECG have a high CC (ranging from 0.65 to 0.95) and lower APL (ranging from 3 to 8) which sets them apart from RNs of chaotic and hyperchaotic systems. We observe that the values of LD are highest for DR and lowest for BB among all ECG. Since the corresponding values of APL are lowest and highest respectively, it suggests that the RNs are more connected for DR and less so for BB. For healthy, high values of CC and APL are observed which may indicate more complexity than unhealthy cases. Moreover, the values of all measures for healthy data are marked by least variability among cases.

4 Scaling of measures with recurrence threshold

As mentioned, the recurrence threshold $\varepsilon$ is crucial in the construction of RNs. In usual practice, it is chosen as the minimum value that makes all the nodes connected as a single network and subsequent analysis is performed on the RN thus constructed. It is clear that increasing $\varepsilon$, will make the network more and more connected. We conjecture that how the network measures change as $\varepsilon$ is varied, can provide a finer level of characterisation, which is not studied so far. We show we can derive relevant information about the nature of the embedded attractor and hence the data, from the nature of variation or scaling of the network measures with increasing $\varepsilon$.

4.1 Variation of degree distribution with recurrence threshold

We vary recurrence threshold from 0.1 (when the networks become just connected) to 1.0, in steps of 0.01. However, we note that there is no significant change after 0.5 and hence we present the results till 0.5 ( Fig.4). In the case of healthy, we observe the bimodal distribution maintains its character but first peak shifts to higher values as we increase $\varepsilon$ while the second peak becomes narrower and increases in height, in addition to shifting to higher degrees. This behaviour tells us that there are at least two spatial scales involved in the underlying dynamics.

In the case of BB, we observe a unimodal distribution at $\varepsilon$ = 0.1, with a single peak (Fig. 3). However, as we increase the threshold to 0.2, a second peak emerges (Fig. 4). For data from CM, we observe a behaviour similar to healthy but the distribution at higher degree saturates earlier. The RNs for data of DR in general show a wider degree distribution. As $\varepsilon$ is increased the second peak becomes narrower but maximum degree is very less as compared to other cases. For RNs from data of MI, the first peak disappears as we increase $\varepsilon$ to 0.5. We have presented the typical behaviour observed for each disease case as distinguished from normal cases. However, because of the presence of secondary diagnoses, some data sets show variability from this typical behaviour reported.

4.2 Scaling of link density with recurrence threshold

The changes in the structure of RN as $\varepsilon$ is varied, can be studied using the link density. This variation as $\varepsilon$ in the range 0.1 to 0.5 is shown in Fig. 5. In general, there exists two scaling regions, corresponding to ranges with two slopes. But the values of slopes and the range of the two scaling regions vary among data sets. The two scaling regions are very pronounced in the case of CM while for MI, there is a gradual change of slope throughout and it is difficult to discern the distinct scaling regions. For BB, on the other hand, the link density is much smaller and there is almost no change of slope. This corresponds to a simple attractor which is very similar to the chaotic Rössler system. In the case of DR the first scaling region is apparent from 0.1 to 0.3, but the second one appears in the range 0.35 to 0.5, corresponding to a sparser attractor with no intermediate structure.

For comparison, we also present results from white noise, Rössler system and Chen system. For Rössler system, we observe an almost linear variation, while for Chen system we observe two scaling regions. However, the change in slope is gradual and the value of link density for smaller thresholds is low, indicating that the regions in the attractor on a small scale are less dense compared to ECG data. The link density for white noise, on the other hand, rises very sharply with $\varepsilon$ and saturates at a very small value of around 0.2, as expected.

The presence of two scaling regions indicates typical structure for underlying attractor of ECG data and the differences in values of slopes can be an indicator to distinguish different cases of disease from healthy ECG. For this, we take the scaling of LD with $\varepsilon$ as LD $\sim$ $\varepsilon^{\gamma}$ and compute values of $\gamma$ from their log-log plots. The results for the indices $\gamma_{1}$ and $\gamma_{2}$ corresponding to the scaling regions are given in Table 2. We would like to emphasise that although the behaviour of ECG data shown in Fig. 4 and 5 are shown for typical cases from each type, the results in tables 1 and 2 are average values computed over all data sets. For standard systems, the values are averaged over data generated from different initial conditions. The presence of secondary diagnosis obscures some results and makes it difficult to separate disease specific features. However, the trends on an average can still be picked up and can guide further, detailed analysis of the time series.

5 Conclusion

In this study, we construct recurrence networks from clinical ECG data of 1-minute duration and analyse their measures for characterisation. We also compare them with those from the standard nonlinear dynamical systems of Rössler (chaotic) and Chen (hyperchaotic), and white noise. The RNs from ECG are found to be different from standard nonlinear dynamical systems and noise in significant ways. The degree distribution of RNs from ECG exhibits a bimodal character and thus form a separate class. We note that this is the first time this novel feature is reported for RNs from time series. We reason that the complex dynamics underlying cardiac system, with structures at two spatial scales for the phase space attractor, is revealed through this unique feature.

RNs from ECG in general have significantly higher value for the clustering coefficient and lower value for the average path length. Among different cases of ECG, Bundle Branch block has values for CC and APL that can be identified separately from all other classes. The healthy show least variability in both measures, and on average have a higher value of both CC and APL.

Our study also indicates that there is relevant information to be extracted from the behaviour of measures as recurrence threshold is varied. In particular, we observe two scaling regions in the link density vs $\varepsilon$ curve for ECG data, that is very different from data of standard chaotic systems and white noise. The scaling indices are calculated for all data sets to quantify these characteristics. The study on how the measures vary as we change the threshold, characterise the dynamics at a finer level. This is the first time such variations are analysed to advantage in the classification of data sets. The analysis presented is to be applied to larger number of data sets so that disease specific ranges of measures and scaling indices can be derived. This work is already in progress and the results will be published elsewhere soon.

The present study closely resembles the tools of multifractal analysis from data in many respects. But in addition to being easier to implement, it has the specific advantage that it requires only short data sets for its implementation. This makes the method faster for practical purposes and more effective for short and nonstationary data and hence can be used for quantifying changes with time for data from the same source.

6 Acknowledgements

One of the authors, Sneha Kachhara acknowledges financial support from Council of Scientific and Industrial Research (CSIR), India. We thank Prof. K.P. Harikrishan for useful discussions.