Investigation of Linear and Nonlinear Properties of a Heartbeat Time Series Using Multiscale Rényi Entropy.

The time series of interbeat intervals of the heart reveals much information about disease and disease progression. An area of intense research has been associated with cardiac autonomic neuropathy (CAN). In this work we have investigated the value of additional information derived from the magnitude, sign and acceleration of the RR intervals. When quantified using an entropy measure, these time series show statistically significant differences between disease classes of Normal, Early CAN and Definite CAN. In addition, pathophysiological characteristics of heartbeat dynamics provide information not only on the change in the system using the first difference but also the magnitude and direction of the change measured by the second difference (acceleration) with respect to sequence length. These additional measures provide disease categories to be discriminated and could prove useful for non-invasive diagnosis and understanding changes in heart rhythm associated with CAN.

1. Introduction

Biological signals, including electrocardiograms (ECG) or the electrical activity of the heart, exhibit complex dynamics which are characterized by nonlinearity and nonstationarity and often include random noise due to movement artefacts [1]. Heartbeat time series associated with health and disease have been extensively investigated, where time and frequency domains, as well as nonlinear methods, are being proposed and summarized in a number of communications [2,3,4,5,6,7,8,9,10,11].

Physiological dynamics of the heartbeat time series change with healthy aging [12,13] and disease, but also during different activities such as sleeping [14,15] and exercise [16,17,18,19,20,21]. Other changes in dynamics can be attributed to pathology including cardiovascular disease and heart failure [22,23,24], diabetes [25], depression [26,27] and Parkinson’s disease [24,28,29]. For all physiological and pathophysiological models of autonomic function, heart rate variability (HRV) is calculated from the cardiac interbeat intervals (IBI) of the time series. All models assume that the extrinsic modulation of the heartbeat by the autonomic nervous system (ANS) and the endocrine system affect HRV by either increasing the interbeat interval (parasympathetic influence), or decreasing the interbeat interval (sympathetic influence), or a combination of both. Disturbance of the ANS modulation by pathophysiological processes, such as oxidative stress, can then lead to the characteristic changes in sympathovagal input to the heart associated with cardiac autonomic neuropathy (CAN).

1.1. Heartbeat Interval Time Series

Non-invasive methods that are independent of patient cooperation are preferable in the diagnosis of CAN, but still require further research to confirm their sensitivity and specificity in stratification of CAN progression. The most common method currently used is heart rate variability analysis [30,31]. HRV is a useful indication of the health of the cardiovascular system, and is commonly used in assessing the regulation of cardiac autonomic function. HRV has been described by a variety of measures such as time domain, frequency domain, and non-linear dynamic (NLD) measures. However, time domain and power spectral density determination are not suitable for the analysis of non-linear and long-range correlated time signals [32]. Application of new signal processing techniques based on NLD, on the other hand, provides supplementary information (i.e., hidden underlying mechanisms) regarding physiological and pathophysiological processes involved in cardiovascular function and pathology. Two components of a time series in particular, being the sign and magnitude, allow further investigation into the characteristics of the time series and have been discussed previously [33,34,35].

1.2. Decomposition of the RR Interval Time Series

The current work is based on Ashkenazy [33], who introduced a decomposition algorithm of the RR interval time series by calculating the beat-to-beat increment or first difference (ΔRR = RRn − RRn−1). The first difference series is then decomposed into the magnitude and sign of the increments (|ΔRRRR| and sign (signΔRR) respectively).

Here, we extend this work by introducing the acceleration, defined as the difference between two successive differences, i.e.,

Velocity is defined as the rate of change in the RR interval length and therefore the first order difference (∆RR). Acceleration is then the second order difference (∆2RR). The second difference or acceleration is a measure of the change of RR points with respect to time and indicates the instantaneous acceleration of the heart rate. We propose that acceleration represents an additional descriptive term for a time series. The scaling properties in sign, magnitude and acceleration can then be analyzed by HRV measures, which define the temporal organization of the original time series. Previously detrended fluctuation analysis (DFA) was applied to the sign and magnitude time series [33]. Here, we analyze, sign, magnitude and acceleration using the multiscale Rényi entropy [36,37,38].

1.3. The Rényi Entropy

Entropy measures can be used to quantify the diversity, uncertainty, or randomness of a system, and are hence considered as beneficial tools for analyzing nonlinear time series, including those of short duration, towards identifying underlying pathology [39,40,41]. Global entropy measures, such as approximate entropy (ApEn) [42] and sample entropy (SampEn) [43], were adapted from the correlation dimension [44,45] and Kolmogorov entropy [46]. RR time series, however, are multifactorial and display multiscale characteristics, and thus neither ApEn nor SampEn are ideal for such types of biosignal processing. Nonlinear, multiscale dynamic systems can however be described by scaling exponents [47], as well as several multiscale measures [48,49]. Rényi entropy has several advantages. The major advantage of Rényi entropy is that it is robust for short time series, nonlinearity and nonstationarity. The Rényi entropy introduced here also has the advantage of addressing how the probabilities are calculated by applying a density method rather than a histogram method, which is the standard for calculation of multiscale entropy [49].

In the current work we use the Rényi entropy, which generalizes the Shannon entropy [50] and is defined as: where p is the probability that a random variable takes a given value out of n values, and α is the order of the entropy measure [50]. H(0) is simply the logarithm of n. As α increases, the measures become more sensitive to the values occurring at higher probability and less to those occurring at lower probability, which provides a picture of the RR length distribution within a signal. The probability, p, can be estimated for any sub-sample of RR intervals, by considering the sub-sample as a point embedded in a multi-dimensional space. The sub-sample is assigned a density measure by evaluating other sub-samples in its vicinity. This addresses the coarse-graining problem for the determination of scaling behavior in biosignal time series inherent in previous applications of the entropy measures by applying a Gaussian kernel [49,51]. The Gaussian kernel is calculated as the sum of all contributions from other RR sub-samples with index j: where σ is the dispersion of the function, and replaces the tolerance as suggested by Costa [48]. We designate the number of RR intervals in the sub-sample as π (not to be confused with the irrational number pi), and use the Euclidean distance measure in π dimensions:

Here, we investigate the efficacy of applying multiscale Rényi entropy as a measure of HRV with respect to the sign, magnitude and the rate of change (acceleration) of the biosignal over time.

2. Methods

2.1. Patient Selection

Heart rate tachograms were obtained from data collected at the Charles Sturt Diabetes Complications Screening Clinic (DiScRi), Australia [52] and were approved by the Charles Sturt University Human Ethics Committee. Written informed consent was obtained from all participants. A 20-min lead II ECG recording was taken from participants attending the clinic, using Powerlab hardware with Chart 7 software (ADInstruments, Sydney) during the morning in an ambient temperature room and after the participants were relaxed. Participants were comparable for age, gender, and heart rate, and after initial screening, those with heart disease, presence of a pacemaker, kidney disease or polypharmacy (including multiple anti-arrhythmic medications) were excluded from the study. The status of CAN was defined using the Cardiac Autonomic Reflex Test battery criteria [53]. Each participant was assigned as either without CAN (71 participants), early CAN (67 participants) or definite CAN (NN participants) [54,55].

2.2. ECG Recording and Obtaining the RR Intervals

From the 20-min RR tachogram, a 10 min segment was selected from the middle in order to remove transient start up artefacts and movement at the end of the recording. The RR intervals were then extracted from this shorter recording, and data were visually verified to not include any missed, extra or misaligned (including ectopic) beat detections. No other information was used in this study. The raw RR interval series for each participant was detrended based on smoothness priors formulation [56]. For the purposes of an initial examination of the RR interval recordings, we have selected one recording from each of these three classes as follows. For each participant class (Normal, Early and Definite), the Standard deviation of RR intervals in the time series were calculated. For each participant, we calculated the difference between this and the Median values of the Standard deviation obtained for all participants of the same class. We selected the RR time series closest to the median for that class. The 10 min RR time series for these representatives are shown in Figure 1. Horizontal scales are the same to allow comparison, but vertical scales are as indicated on each graph. Figure 1 shows that the participant from the Normal class manifested RR intervals with mostly low deviation from the mean, but some large excursions (standard deviation = 0.0357). In comparison, the participant from the Definite class showed fewer large excursions (SD = 0.0173), while the participant from the Early class was in between these (SD = 0.025926).

Figure 1

RR interval time series of normal, early cardiac autonomy neuropathy (eCAN) and definite CAN (dCAN).

2.3. Decomposition

The RR interval time series was decomposed into increment, magnitude, sign and acceleration, as discussed above. Raw RR intervals were filtered and the trend was removed. Increments were calculated as the difference between successive RR intervals. The magnitude, sign and acceleration of the increments were then determined. Finally, the Rényi entropy was calculated for the sign, magnitude and acceleration time series, using a variety of values for parameter sequence length π, exponent α, and width of the kernel function σ. This results in four different measures:

Rényi entropy calculated from a sequence of the magnitude of the difference in RR intervals

Rényi entropy calculated from a sequence of the sign of the difference in RR intervals

Rényi entropy calculated from a sequence of the acceleration of RR intervals

2.4. Calculating the Multiscale Rényi (MSRen) Entropy

The Rényi entropy was calculated for scaling exponents α of integer values from −5 to +5. The entropy values were then normalized by dividing by log2 of the number of length of the RR interval time series. A range of sequence lengths, π, was also used, and the dispersion of the Gaussian function (σ) was varied in proportion. Sequence lengths of 1, 2, 4, 8 and 16 RR intervals were adopted, with corresponding values of σ as 0.01, 0.02, 0.04, 0.08 and 0.16, respectively. A Mann-Whitney test was performed to compare the Rényi value obtained for the Normal to that obtained for the Early CAN group, and a similar comparison between the Early CAN and the Definite CAN group, and between the Definite CAN and the Normal group.

3. Results

For each patient group, we calculated the median value of the standard deviation of the RR intervals and selected the patients with standard deviation of RR intervals closest to the median for the group. The resulting three representative patients are used to illustrate the differences in sign, magnitude and acceleration between Normal, Early CAN and Definite CAN. Figure 2, Figure 3 and Figure 4 present a sample of 100 filtered RR intervals and their decomposition, using data from the three representative groups to illustrate the effect of working with the sign of the RR interval, first and second difference. All vertical axes are numbered in seconds. It can be observed, for example, that the increment ΔRR (b) varies between ±0.2 s with excursions up to 0.5 s, while the acceleration (e) varies between ±0.7 ms. There are frequent reversals of sign (d) with some periods of a continuation of the same sign.

Figure 2

An illustration of the composition of the raw 100 RR tachogram for participants without CAN (Normal group). (a) The RR time series after filtering and pre-processing; (b) The increment (ΔRR = RRn − RRn−1) of the time series shown in (a); (c) The magnitude of the increment; (d) The sign of the increment; (e) The acceleration (Δ2RR = RRn − 2RRn−1 + RRn−2).

Figure 3

An illustration of the composition of the raw 100 RR tachogram for a participant with Early CAN. (a) The RR time series after filtering and pre-processing; (b) The increment (ΔRR = RRn − RRn−1) of the time series shown in (a); (c) The magnitude of the increment; (d) The sign of the increment; (e) The acceleration (Δ2RR = RRn − 2RRn−1 + RRn−2).

Figure 4

An illustration of the composition of the raw 100 RR tachogram for a person with Definite CAN. (a) The RR time series after filtering and pre-processing; (b) The increment (ΔRR = RRn − RRn−1) of the time series shown in (a); (c) The magnitude of the increment; (d) The sign of the increment; (e) The acceleration (Δ2RR = RRn − 2RRn−1 + RRn−2).

Figure 3 shows similar information from the representative participant with early CAN. The range of variation in RR interval, ΔRR and |ΔRR| can be observed to be much smaller than those in Figure 2, indicating a smaller variance in the RR interval. In addition, the acceleration is large compared to the representative with early CAN in Figure 3. The sign (Figure 2e) shows fewer changes in direction compared to the representative with early CAN in Figure 3.

Figure 4 shows a sample of the information from a participant with definite CAN. The difference in RR intervals (Figure 4b,c) can be seen to be even smaller than the example shown in Figure 2 or Figure 3, while the acceleration (Figure 4e) also has a smaller range than the example from either the Normal or Early group. The sign is different to either of the previous examples, as there appears to be more frequent reversals in sign when compared to those examples, but there is less of a mixture of fast and slow changes in sign.

In order to quantify these differences, the variation of each time series was evaluated, for each participant in the study, using the Rényi entropy. A variety of values were used for the parameters (sequence length π, exponent α and width of the kernel function σ).

Our results comparing Normal (N), Early (E) and Definite (D) CAN, based on the magnitude of the difference of RR intervals for sequence length π set to 1, 2, 4, 8, 16, and for α = +5 applying the Mann-Whitney tests obtained the smallest p-value (p < 0.0001) for the Definite to Normal comparison with a sequence length of π = 1 values for ΔRR. Normal versus Early was best differentiated with longer sequences of length π = 4 or π = 8, whereas for Early versus Definite the optimal sequence length was again π = 1. Extracting the sign of a ΔRR sequence provides information on the linear aspects of the traces but the separation of the classes is less pronounced as seen in the figures above and results in the tables below and hence the p-values are much larger, indicating a lesser role of the linear characteristics of the signals in differentiating CAN progression. Only Normal to Early CAN was significantly different for a sequence length of π = 8, suggesting that the nonlinear, fractal-like characteristics may play a larger role in CAN development.

Separating the classes based on the acceleration of ΔRR increments results in the smallest p-value (8.13 × 10−5) obtained for the Mann-Whitney test comparing Definite CAN to Normal, and a sequence length π = 4. For acceleration, separation of CAN progression improves with sequence length up to π = 4, and then decreases again for all comparisons, except for Normal versus Early, where the best separation is seen using π = 16. However, the best overall comparative results were found with π = 4.

Table 1 concerns the magnitude of differences (|ΔRR|) and shows p-values obtained from three Mann-Whitney tests comparing Normal to Early (NE), Early to Definite (ED) and Definite to Normal (DN), for different values of the Rényi parameters sequence length π, and exponent α. The width of the kernel function σ was always chosen as π/100. Nearly all of these p-values are significant at the p < 0.01 level, and some are extremely small, suggesting that these Rényi exponents show an effect for all three of these comparisons.

Table 1

Classification results based on Rényi exponents applied to the magnitude of differences (|ΔRR|). Figures represent p-values for the results of Mann-Whitney tests for comparisons of Normal to Early (NE), Early to Definite (ED) and Definite to Normal (DN). Values shown in bold are the smallest p-value for each comparison (Table 1), while tests that were not significant are indicated by n.s.

	Test	π = 1	π = 2	π = 4	π = 8	π = 16
α = 1	NE	0.0001	<0.0001	<0.0001	0.0001	0.002
	ED	0.004	0.006	0.02	n.s.	n.s.
	DN	<0.0001	<0.0001	0.0002	0.002	0.03
α = 2	NE	0.0002	<0.0001	<0.0001	<0.0001	0.0007
	ED	0.002	0.004	0.01	0.05	n.s.
	DN	<0.0001	<0.0001	0.0001	0.001	0.02
α = 3	NE	0.0002	<0.0001	<0.0001	<0.0001	0.0004
	ED	0.003	0.005	0.009	0.004	0.2
	DN	<0.0001	<0.0001	<0.0001	0.001	0.01
α = 4	NE	0.0003	<0.0001	<0.0001	<0.0001	0.0003
	ED	0.002	0.005	0.009	n.s.	n.s.
	DN	<0.0001	<0.0001	<0.0001	0.0009	0.007
α = 5	NE	0.0003	<0.0001	<0.0001	<0.0001	0.0002
	ED	0.002	0.005	0.009	0.02	n.s.
	DN	<0.0001	<0.0001	<0.0001	0.0007	0.006

In general, the smallest values are found for normal versus definite CAN as would be expected. However, the table suggests that some values of the Rényi parameters are better than others at demonstrating this effect. Generally, the sequence length of 2 provides the best separation for ED (early–definite) and DN (definite–normal), but using π = 4 provides the best separation for NE (normal–early). n.s.—not significant.

Table 2 illustrates results for acceleration and shows p-values obtained from three Mann-Whitney tests for different values of sequence length π, and exponent α. The figures show an optimum value for π = 4 and for a variety of values for α. For short sequence lengths, the p-value increases with increasing α. For longer sequences the opposite is true.

Table 2

Classification results based on Rényi exponents applied to acceleration (|Δ2RR|). Figures represent p-values for the results of Mann-Whitney tests for comparisons of Normal to Early (NE), Early to Definite (ED) and Definite to Normal (DN). Values shown in bold are the smallest p-value for each comparison (Table 1), while tests that were not significant are indicated by n.s.

	Test	π = 1	π = 2	π = 4	π = 8	π = 16
α = 1	NE	0.003	0.0005	0.0002	0.0003	0.0005
	ED	n.s.	0.03	0.01	0.02	n.s.
	DN	0.00147	0.0002	<0.0001	0.0005	0.004
α = 2	NE	0.007	0.0008	0.0003	0.0004	0.0004
	ED	n.s.	0.02	0.009	0.02	0.05
	DN	0.002	0.0003	<0.0001	0.0003	0.001
α = 3	NE	0.01	0.001	0.0003	0.0005	0.0003
	ED	n.s.	0.03	0.007	0.02	0.04
	DN	0.002	0.0004	<0.0001	0.0002	0.001
α = 4	NE	0.01	0.0009	0.0004	0.0005	0.0003
	ED	n.s.	0.03	0.007	0.02	0.03
	DN	0.003	0.0004	<0.0001	0.0002	0.0009
α = 5	NE	0.01	0.001	0.0004	0.0006	0.0002
	ED	n.s.	0.03	0.007	0.02	0.03
	DN	0.0038	0.0004	<0.0001	0.0001	0.0006

The actual values of the Rényi entropy calculated from the magnitude of the increment of RR intervals |ΔRR|, using the parameters sequence length π = 4, and width of the kernel function σ = 0.04, are illustrated in Figure 5. The exponent α was varied so that −5 ≤ α ≤ 5. The inset highlights details of the exponents corresponding to the positive values of α. Rényi entropy calculated for the class of Early CAN lies in between those for the Normal and Definite classes.

Figure 5

Values of Rényi entropy based on the magnitude of the difference between RR intervals, for sequences of 4 values.

Rényi entropy calculated from the acceleration of the RR intervals Δ2RR, using the parameters sequence length π = 4, and width of the kernel function σ = 0.04 indicates a better separation for positive α (Figure 6.).

Figure 6

Values of Rényi entropy based on the acceleration of RR intervals, for sequences of 4 values.

4. Discussion and Conclusions

In physiological dynamic systems, various extended concepts of entropy, such as approximate entropy (ApEn), sample entropy (SampEn), and multi-scale entropy have been developed to quantify various physiological signals [42,43,57]. The major advantage of Rényi entropy is that it is robust for short time series, nonlinearity and nonstationarity. The Rényi entropy introduced here also has the advantage of addressing how the probabilities are calculated by applying a density method based on a Gaussian kernel rather than a histogram method, which is the standard for calculation of multiscale entropy. H(α) is the order of the Rényi entropy measure. As α increases, the measures become more sensitive to the values occurring at higher probability and less to those occurring at lower probability, which provides a picture of the RR length distribution within a signal [49]. Including acceleration in our model then adds information about the heart rate variability by providing information not only on the change in the system using the first difference but also the magnitude and direction of the change measured by the second difference (acceleration) with respect to sequence length.

Heart rhythm is characterized by a scale-invariant, nonlinear dynamics displaying long-range power-law correlations over a range of time scales [58,59] akin to 1/f [60] or fractal-like scaling [1,39,47,57,58,61,62,63,64,65]. Fractal-like scaling analysis has been shown to indicate risk of adverse cardiac events [66,67].

Bio signals quantifying cardiac interbeat intervals (RR intervals) exhibit complex dynamics that vary with age and disease and can be characterized by scaling laws [12,68]. In healthy subjects, RR interval time signals present large variability, which is a function of the numerous physiological processes that influence heart rhythm, including ANS and neuroendocrine factors. Nonstationarity and nonlinear dynamics characteristic of these signals are believed to be due to the complex interaction between the two branches of the ANS, endocrine factors and the intrinsic cardiac control mechanisms.

Early identification of CAN is crucial for more effective clinical outcomes. Studies have shown that the one of the earliest signs for a CAN diagnosis is the reduction of HRV. Thus, understanding of the time series characteristics and selecting an appropriate method to analyze these signals and interpret the results is paramount. A consistent finding of ours is that the most difficult two classes to separate were Definite and Early CAN. This implies that patients in the early stages of CAN have similar HRV features to those in the definite group. This may be a reflection that the existing CART criteria are somewhat conservative in identifying CAN, or the two blood pressure tests included in the CART battery indicative of sympathetic dysfunction do not clearly identify disease progression from early to definite CAN, and that sympathetic dysfunction may already be a factor in early CAN [69].

The typical RR tachogram consists of linear and nonlinear portions, which overlap and lead to the characteristic heart rate variability. In this work, we show that different components of the RR tachogram are able to differentiate between the stages of CAN progression from normal and early CAN to definite CAN. These different components rely on the fact that control of heart rate entails changes in both a positive and negative directions. In particular, the magnitude and acceleration of the changes in RR increments separate all three groups. Both of these series carry information on the nonlinear properties of the interbeat interval time series and indicate that fractal-like or power law dynamics within the biosignals become more prominent with disease progression. This complex behavior is further illustrated by the larger and more often occurring deviations in acceleration.

Recent NLD methods continue to shed light on HRV changes under various physiological and pathological conditions, providing valuable potential prognostic and diagnostic information and complementing traditional time- and frequency-domain analyses. With the advent of multiple tools and algorithms, it is critical to identify which of these methods should be selected and under which conditions they should be applied. Our work aligns with previous work, confirming the efficacy of complex measures for representing and quantifying heart rate variability. The current research has focused on investigating the differences in successive RR intervals adopting interbeat acceleration as a novel feature, which provides additional information about the nonlinearity of heartbeat regulation and hence the identification of disease. The high degree of separation obtained between classes of disease points to its diagnostic and risk stratification potential in cardiac autonomic neuropathy, and provides a much less invasive test for this disease, with the advantages of faster diagnosis, better access to treatment and more effective clinical outcomes.