1 Introduction

Cardiovascular autonomic neuropathy (CAN) is highly prevalent and a serious complication in patients with diabetes mellitus. CAN progression may lead to severe postural hypotension, exercise intolerance, enhanced intraoperative instability, increased incidence of silent myocardial infarction and ischemia [36]. Around 75 % of people with diabetes die from cardiovascular disease such as heart attack and stroke, which includes autonomic neuropathy as a cause [17, 33]. Early detection of CAN in diabetic patients and intervention are therefore of prime importance to reduce the increased mortality of diabetes patients.

A non-invasive autonomic nervous system test battery, consisting of five tests, was specifically designed for identifying CAN by Ewing [9]. The major drawbacks of the Ewing battery are: (1) it requires patient cooperation [9]; (2) it is unable to be performed when co-morbidities are present in the patients such as heart or lung disease [25]; (3) it is less sensitive to changes associated with cardiac autonomic neuropathy compared to spectral methods [8].

Since the change in heart rate variability (HRV) is regarded as one of the early signs of CAN [32], the most common method used for CAN diagnosis is HRV analysis [29]. However, conventionally used time and frequency domain parameters of HRV are not always suitable for analysis because of the non-stationarity characteristic of the ECG recordings and the presence of nonlinear phenomena in the physiological signal’s parameter variability [10]. Spectral methods also lack sensitivity, especially for detecting asymptomatic/preclinical CAN [14]. Therefore, new methods that are non-invasive and independent of patient cooperation are preferable in the diagnosis of preclinical CAN. Only a few studies have applied new parameters based on nonlinear dynamics theory to HRV analysis in diabetic mellitus (DM) patients [7, 10, 13, 15, 37]. However, the maximum ROC area reported by Khandoker et al. was 0.79 when using sample entropy (SampEn) for a study of 17 subjects (9 CAN+ and 8 CAN−). The association of CAN with the alteration of sympatho-vagal balance, reported by Khandoker et al. [16], showed that tone–entropy (T–E) could be a better marker for identifying CAN. An important benefit of using T–E analysis is that it is not influenced by the time period of data acquisition, nor the baseline heart rate. In addition, T–E data processing has no signal deformation process such as filtering, window size based on window lengths and there is no need to control respiration rate, allowing data to be analyzed as it is recorded [23].

The conventional T–E method of quantifying HRV uses successive RR intervals with the implicit assumption that the current beat is influenced by the immediately preceding beat. However, the delay of such influence is affected by the variation of baroreflex sensitivity and therefore can be greater than one beat [4, 24]. For example, the baroreflex sensitivity is reduced in chronic renal failure (CRF) patients [34], which augments the delay in the heart rate response. Therefore, a heart beat influences not only the beat immediately following it, but also up to 6–10 beats downstream [19]. Lerma was the first to show this with several researchers confirming higher lag being associated with disease processes [19, 21]. Thus, successive RR interval duplets will likely underestimate the role of the autocovariance function of RR intervals i.e., the ability of heart beats to influence a train of succeeding beats. Moreover, the autocovariance function of RR intervals captures additional aspects of HRV (e.g., non-linearity) that are otherwise masked by the strong correlation between successive beats if lag 1 (n vs. n + 1 beats) tone–entropy is used. Therefore, we hypothesized that multi-lag T–E analysis can overcome the limitations of the present practice of single lag T–E analysis in HRV studies.

The aims of this study were (1) to study the impact of data length on T–E values; (2) to study the T–E process at multiple lags; (3) to define parameter values (length and lag) for calculating multi-lag T–E of HRV signals; and (4) to compare the performance of multi-lag T–E with conventional T–E measurement in identifying CAN using HRV signal during supine rest.

2 Methods

2.1 Multi-lag T–E analysis of HRV signal

A RR interval or period is defined as the time difference between two consecutive R peaks of the electrocardiogram (ECG) signal. Let the RR intervals time series RR be defined as:

$$ {\text{RR}} \equiv \left( {{\text{RR}}_{1} ,\;{\text{RR}}_{2} , \ldots ,{\text{RR}}_{N} } \right) $$

where, N is the number of RR intervals. Heart rate acceleration and inhibition can be determined from the difference of consecutive RR intervals. If RRi+1 become shorter than RRi then it is an acceleration of heart rate. Therefore, acceleration of the heart is expressed as a plus difference and inhibition as a minus difference of RR intervals. However, to reduce the impact of heart rate variation over a wide range of time and different subjects, normalized variation in RR interval is preferred to monitor the variability. In conventional T–E analysis, percentile change of the successive RR intervals with respect to the first RR interval is expressed as the percentage index (PI) and defined as:

$$ {\text{PI}}\left( i \right) = \frac{{{\text{RR}}_{i} - {\text{RR}}_{i + 1} }}{{{\text{RR}}_{i} }} \times 100 $$
(1)

The tone is defined as a first order moment (arithmetic average) of this PI time series as:

$$ {\text{Tone}} = \frac{1}{N - 1}\sum\limits_{i = 1}^{N - 1} {{\text{PI}}\left( i \right)} $$
(2)

Tone is the balance between accelerations (PI > 0) and inhibitions (PI < 0) of the heart rate and represents the sympatho-vagal balance faithfully as shown in previous studies [2, 22]. Entropy is defined from the probability distribution of PI using Shannon’s formula [30]:

$$ {\text{Entropy}} = - \sum\limits_{i = 1}^{n} {p\left( i \right)\log_{2} p\left( i \right)} $$
(3)

where, p(i) is a probability of PI having values in the range i < PI < i + 1, where i is an integer. The entropy evaluates total acceleration–inhibition activities, or total heart period variations, in a familiar unit of bit.

For multi-lag T–E analysis, we have introduced the lag (m) in Eq. (1), used to derive the PI time series from the RR time series signal. Hence, in the multi-lag T–E analysis PI is expressed as the percentile change of the i-th and i + m-th RR intervals with respect to the i-th RR interval and is defined as:

$$ {\text{PI}}\left( i \right) = \frac{{{\text{RR}}_{i} - {\text{RR}}_{i + m} }}{{{\text{RR}}_{i} }} \times 100 $$
(4)

where, m is an integer and m = 1 represents the conventional T–E analysis. The detailed methodology of conventional T–E analysis has been described in previous reports [2, 22].

An important benefit of conventional (lag-1) T–E evaluation is that it is not influenced by the time period of data acquisition and increasing lag time reintroduces the sensitivity to recording length. Therefore, we analyzed the variation of tone and entropy values with varied length of the RR interval data (from 50 to 900 beats) and a range of lags (1 ≤ m ≤ 8).

2.2 Patients and ECG signals

After standard exclusion criteria were applied to ensure that any changes in HRV detected were due to the severity of cardiac autonomic neuropathy, 41 patients with Type 2 diabetes mellitus were included in the study. Fifteen patients had definite CAN (CAN+), whilst the remaining 26 were negative for CAN (CAN−), being without clinical signs and symptoms of CAN. Detailed demography of subjects is shown in Table 1.

Table 1 Subject demography of CAN+ and CAN− groups
Full size table

The research protocol was approved by the Charles Sturt University Ethics in Human Research Committee (03/164) and complies with the declaration of Helsinki. All patients gave written informed consent for this study. Exclusion criteria for the CAN− group included those with a history of cardiac pathology, hypertension or use of antihypertensive or antiarrhythmic medication and those with less than 85 % qualified sinus beats. For the CAN+ group all participants with greater than 85 % qualified sinus beats were included. ECGs were recorded over 20 min using a lead II configuration (Maclab ADInstruments, Australia) and recorded on Macintosh Chart Version 5 with a sampling rate set at 400 Hz and a notch filter at 50 Hz. ECG signals were edited using the MLS310 HRV module (version 1.0, ADInstruments, Australia) included with the Chart software package. High frequency noise was removed with a 45 Hz low pass filter and a 3 Hz high pass filter adjusted for wandering baseline. RR intervals, and the difference of successive R waves of the QRS complex, were calculated using the algorithm developed by Pan and Tomkin [26]. Ectopic beats were selected visually and deleted manually. CAN+ was determined using the complete Ewing battery [9].

2.3 ROC area analysis

In order to retain the relative importance of features, receiver-operating curve (ROC) analysis was used [11], with the area under the curve for each feature represented by the ROC area. A ROC area value of 0.5 indicates that the distributions of the features are similar in the two groups with no discriminatory power. Conversely, a ROC area value of 1.0 would mean that the distribution of the features of the two groups do not overlap at all. The ROC is obtained by automatic selection of different thresholds or cut-off points and calculating the sensitivity/specificity pair for each one of the cut-off points using MATLAB. The area under the ROC curve was approximated numerically using the trapezoidal rules [11] where the larger the ROC area is, the better the discriminatory performance.

2.4 Classification of CAN patients

In this study, we used the quadratic discriminant (QD) classifier to test the ability of tone and entropy values together in detecting CAN+ subjects. The beat sequence length len = 250 and lag m = {1, 2, 3} was taken for the tone and entropy calculation based on the ROC areas. A leave-one-out (LOO) cross-validation scheme was adopted to evaluate the generalization ability of the classifiers. Cross-validation procedures have been used in a number of classification evaluations, particularly for limited data sets [28]. In this scheme, the data set was uniformly divided into 41 subsets with one used for testing (unknown CAN classification) and the remaining 40 records used to train the classifiers. This was repeated for the remaining subsets so that all subsets were used as the testing sample.

The following three measures of accuracy, sensitivity, and specificity were used to assess the performance of the classifiers [5, 27]:

$$ \begin{array}{lll} {{\text{Accuracy}} = \frac{{{\text{TP}} + {\text{TN}}}}{{{\text{TP}} + {\text{FP}} + {\text{TN}} + {\text{FN}}}} \times 100} \\ {{\text{Sensitivity}} = \frac{\text{TP}}{{{\text{TP}} + {\text{FN}}}} \times 100} \\ {{\text{Specificity}} = \frac{\text{TN}}{{{\text{TN}} + {\text{FP}}}} \times 100} \\ \end{array} $$

where, TP is the number of true positives, i.e., the classifier identifies a patient that was labeled as CAN+; TN is the number of true negatives, i.e., the classifier identifies a patient that was labeled as CAN−; FP is false CAN+ identifications; and FN is false CAN− identifications. Accuracy indicates overall detection accuracy, sensitivity is defined as the ability of the classifier to accurately recognize a CAN+, whereas specificity indicates the classifier’s ability not to generate a false negative (CAN−).

2.5 Statistics

In this study, we used an ANOVA analysis assuming unknown and different variance for testing the hypothesis regarding mean i.e., the mean of tone values for data length len = 900 and other lengths len = 50–850 groups are equal. The same test has been used to test the hypothesis for entropy values.

Non-parametric Kruskal–Wallis has been used to investigate the significance between age distributions of subjects of CAN− and CAN+ group. Since, the number of subjects are small and their distribution is not normal non-parametric test is more appropriate than parametric test.

The effect size—in essence how much overlap is there between two groups—was investigated using Cohen’s d value [6] and then the power analysis, which is the probability that the null hypothesis (there is no difference in T–E values between CAN− and CAN+ groups) is wrong, is carried out to justify the repeatability of the difference in T–E values found between two groups for small sample size. Large Cohen’s d value represents less degree of overlap between two groups and vice versa. On the other hand, a small p value indicates that the null hypothesis can be rejected.

3 Results

Mean tone and entropy values for both the CAN− and CAN+ group for all lags (1 ≤ m ≤ 8) and beat sequence lengths (50 ≤ len ≤ 900) are shown in Fig. 1. Mean tone values are lower for all lags and beat sequence lengths in the CAN− compared to the CAN+ group. In addition, mean tone values in the CAN− group consistently decrease with increasing lag for all beat sequence lengths used in this study. In contrast, the mean tone values associated with a specific lag tend to increase with increasing length of the ECG analyzed (not consistent at every lag). Similarly, the mean tone values in the CAN+ group decrease with increasing lag for any sequence length greater than 150 beats. However, there was no consistent relationship between the mean tone values at a specific lag and sequence length.

Fig. 1
figure 1

Mean tone and entropy values for the CAN− and CAN+ group for all lags (1 ≤ m ≤ 8) and beat sequence lengths (50 ≤ len ≤ 900)

Full size image

Mean entropy values are higher in the CAN− than CAN+ group for all lags and beat sequence lengths (Fig. 1). The mean entropy values in the CAN− group are mostly increased (except m = {6, 7} and len = 50; m = 6 and len = 500) with increasing lags. Similarly, the mean entropy values for a specific lag mostly increased with increasing recording length of the ECG (beat sequence length). The mean entropy values for the CAN+ group increased (except lag = {8} and len = 50) with increasing lags, whilst the mean entropy value for a specific lag showed an increasing trend only up to len = 200 where this trend disappeared.

The association between T–E related estimates of HRV and beat sequence lengths for varying lags was analyzed and shown in Fig. 2 for the associations of two extreme examples of beat sequence lengths (len = 50, Fig. 2a, b; len = 900, Fig. 2c, d). Similarly, we analyzed the association between T–E related estimates of HRV and lag for varying beat sequence lengths. Figure 3 shows these associations for two lags.

Fig. 2
figure 2

Tone–entropy values in the CAN− and CAN+ group associated with beat sequence length. Panels a–b show tone–entropy with respect to increasing length at lag-1 and panels c–d show tone–entropy for lag-8

Full size image
Fig. 3
figure 3

Tone–entropy response to changing lag in patients with definite CAN and negative CAN. Panels a–b show the T–E response for sequences 50 beats long, whereas panels c–d show the T–E response for sequences 900 beats long

Full size image

ROC areas are calculated to measure the performance of tone and entropy features in differentiating between CAN− and CAN+ for all lags and beat sequence lengths. ROC areas of tone for all lags and beat sequence lengths len ≥ 250 are summarized in Table 2. The maximum ROC area (=0.95) between CAN− and CAN+ using tone is found for m = 2 and len = {450 and 800}. For any beat sequence lengths, the maximum ROC area for classifying the CAN− and CAN+ groups is found for lag m = 2. Beat sequence length of len = 400 for lag m = 3 also resulted in a maximum ROC.

Table 2 ROC areas of tone and entropy for different lags and beat sequence lengths
Full size table

ROC areas of entropy for all lags and beat sequence lengths len ≥ 250 are summarized in Table 2. The hypothesis and method for selecting length len ≥ 250 are discussed in detail in Sect. 4.2. The maximum ROC area (=0.97) between CAN− and CAN+ using entropy is found for m = 3, len = 750 and m = {2, 3} and len = {800, 850 and 900}. Maximal ROC results were also found for other combinations of length and lag as shown in Table 2.

In this study, we have used the quadratic discriminant with LOO testing methodology to quantify the accuracy of tone and entropy values in discriminating CAN+ from CAN−. LOO allows determination of how accurate an unknown recording is classified into the correct class (CAN− or CAN+). The beat sequence length len = 250 was used for classification, as T–E values become consistent over multiple lags at length len ≥ 250 (see Sect. 4.2). In addition, lag m = 1–3 was selected as the maximum ROC area between CAN− and CAN+ was found at lag 2 and 3 (see Table 2) and the result needs to be compared to lag 1 or conventional T–E analysis. Results of the LOO cross-validation tests (accuracy, sensitivity and specificity) of the QD classifier are summarized in Fig. 4 and Table 3. For beat sequence length len = 250, the accuracy reached 100 % at a lag m = {2, 3}. The effect size and power (one-tailed and two-tailed) for tone and entropy parameters at lag m = {1–3} and for length len = 250 are shown in Table 4.

Fig. 4
figure 4

Classification performance of the quadratic classifier between CAN− and CAN+ subjects using tone and entropy features with beat sequence lengths len = 250 and lags m = {1, 2, 3}

Full size image
Table 3 Sensitivity, specificity and accuracy obtained using QD classifier with tone and entropy feature to classify CAN+ from CAN− using the minimum beat sequence length len = 250 and lags m = {1, 2, 3}
Full size table
Table 4 Effect size (Cohen’s d) and power (one-tailed and two-tailed) values of tone and entropy features of CAN− and CAN+ subjects for beat sequence length len = 250 and lags m = {1, 2, 3}
Full size table

4 Discussion

Linear and nonlinear analysis of RR intervals has been applied for some time for the classification of cardiovascular disease and cardiac autonomic neuropathy in particular. These methods differ in terms of their appropriateness for analysis of the RR interval time series, which is non-stationary and nonlinear and also with respect to the information they provide including their sensitivity and specificity [1, 7, 8, 12, 18, 20]. For many of the current methods, the assumption is that for all consecutive beats, each beat only has influence on the subsequent beat, whereas it has been shown that each beat in a time series can have an influence on up to ten beats downstream [19]. This opened up the question of finding a method that allows measuring this multi-lag characteristic of RR interval time series and its usefulness in clinical application.

Our paper confirmed the finding of [19] that a current beat can influence a beat further downstream and extended this work by applying multi-lag analysis to cardiac autonomic neuropathy. Our approach using multi-lag T–E values of HRV can be used to correctly identify a patient as CAN+ or CAN−, which takes into consideration the influence an individual beat has on subsequent beats in the time series. The physiological interpretations of tone and entropy in various experimental settings were previously reported [16, 22]. Lower tone values (negative) indicate that vagal activity predominates in the sympatho-vagal balance in a healthy population at rest [3, 22, 35]. This is also reflected in higher entropy values in the healthy (CAN−) population. Higher tone (nearing zero) and lower entropy values indicate that parasympathetic efferent pathways are progressively withdrawing their activity [16]. The major findings of the study are discussed in the following subsections.

4.1 Changes of T–E values at different lags

For any beat sequence length (50 ≤ len ≤ 900), the mean tone value for the CAN− group decreased with increasing lag. This indicates that the parasympathetic influence of the current beat decreases with increasing lag (distant heart beats) in the CAN− group. A similar response was found for the CAN+ group with a minimum beat sequence length of 200 beats. Therefore, tone can be considered as a measurement of sympatho-vagal balance at multiple lags with the early lag (~1–5 beats) predominantly indicating parasympathetic influence.

For any beat sequence length (50 ≤ len ≤ 900), the mean entropy values in the CAN− group increased (except m = {6, 7} and len = 50; m = 6, and len = 500) with increasing lag (Fig. 1). Higher entropy values indicate that the degree of both sympathetic and parasympathetic activity increases with increasing lag in CAN− patients. Similarly, the mean entropy values in the CAN+ group tend to increase (except lag = {8} and len = 50) with increasing lag. The increase in entropy with increasing lag for both the CAN− and CAN+ group suggests that the correlation between parasympathetic and sympathetic function decreases with increasing lag. Our results show that the parasympathetic influence has a range up to lag 5 and then decreases as the sympathetic influence increases. Therefore, at higher lag the variance of entropy decreases and becomes more stable as it mainly reflects the sympathetic and endocrine components.

The influence of a heart beat on heart beats downstream from it may be a function of sinus arrhythmia as proposed by [19, 31] with Poincaré Plot analysis. Using Poincaré analysis, HRV measures, which are influenced by respiratory sinus arrhythmia, change their pattern at a lag value of 6 or 7 [31]. However, since in our study the pattern of tone or entropy values did not change (from increasing to decreasing or vice versa), physiological factors other than respiratory sinus arrhythmia may influence T–E values at higher lags. Alternatively, this finding may be a function of the different HRV analysis used here compared to the Poincaré method described in [19]. This requires further investigation by evaluating the baroreflex activity in conjunction with multi-lag T–E analysis and comparing the results with different HRV algorithms.

4.2 Changes of T–E values with incremental beat sequence lengths

We have observed the association between the T–E estimation and varying the length of beat sequences. Conventional T–E analysis (lag-1) is reported as an evaluation process, which is not influenced by the time period of data acquisition. We have demonstrated here the effect of varying beat sequence length on T–E values at higher lags when using the multi-lag T–E algorithm.

The distribution of both tone and entropy values changes with length for both CAN− and CAN+ groups. This change can be either due to the length/number of beats used for analysis or due to changes in function of the autonomic nervous system. To reduce the influence of beat sequence length on the tone and entropy parameters, we have defined the minimum length for reliable T–E analysis. The hypothesis behind the minimum beat sequence length is that the minimum length should be the length for which the mean of the distribution of the parameters (tone and entropy) is insignificantly different to the mean of the distribution for maximum data length. We have tested the hypothesis using an analysis of variance (ANOVA) between T–E values at beat length 900 (maximum for this study) and beat lengths down to 50 beats. The ANOVA analysis was performed for all lag. For the CAN− group, the tone value is insignificantly different for beat lengths ≥250 compared to shorter lengths. In contrast, the entropy values in the CAN+ group at length 50 showed the same mean distribution with length 900. Therefore, we selected a minimum length of 250 beats to make the T–E measurement consistent.

5 Detection of CAN in diabetic subjects using T–E analysis

It should be noted that the mean tone values are lower in CAN− than in the CAN+ group for all lags and beat sequence lengths. From a pathophysiological perspective this indicates that the parasympathetic predominance found in normal beat regulation is reduced in CAN+ patients to some extent. Our finding supports our previous results that heart function with respect to CAN+ resembled that of the parasympathetically denervated heart [16]. The change in the predominance of the parasympathetic activity, measured using tone, is best captured at lag m = 2 or m = 3 rather than m = 1.

Another important finding is that the mean entropy values are higher in the CAN− group than in the CAN+ group for all lags and beat sequence length. Entropy is a measurement of total autonomic activity [16, 22] and in CAN+ subjects this total activity consisting of the sympathetic and parasympathetic influence is less than in CAN− subjects. However, the reduction in entropy values in the CAN+ group may be attributed to either a withdrawal of vagal activity, which is reflected by the increased tone or increased sympathetic activity. This is again consistent with our previous findings [16, 22, 31, 35]. From a pathophysiologial perspective the higher activity of the parasympathetic system in maintaining a steady heart rate and the associated higher firing rate of the parasympathetic nerves increases the likelihood of free radical damage of the nerve terminals and a reduction in the magnitude of the parasympathetic influence on heart rate. The change in neural activity measured by entropy is best captured at lag m = {2, 3} leading to higher accuracy when classifying CAN using entropy.

Use of multi-lag T–E analysis introduces a novel analysis of cardiovascular function with respect to autonomic modulation and can be used to obtain better results in classifying CAN− and CAN+ patients. From the results of this study, we can conclude that CAN− and CAN+ can be better differentiated at lag m = {2, 3} than at lag m = 1. Moreover, the data length used for multi-lag T–E analysis is 250 beats (ECG signal of less than 5 min), which is shorter than the minimum data length required of many time and frequency domain analysis. CAN manifests as a deterioration of autonomic modulation at 2 or 3 beats, which is in line with the findings that the parasympathetic response of heart rate occurs within 2–5 beats. The results of this study suggest that multi-lag T–E is a sensitive indicator of sympatho-vagal balance and activity that may be helpful in detecting CAN.

Although the number of subjects are small in each group of this study, Cohen’s d value was large for length len = 250 and lags m = {1–3}, which indicates that the two groups are substantially different and is reflected in the power (Table 4). Therefore, we conclude that the difference in T–E values among CAN− and CAN+ groups, as well as the classification results found in this study can be repeated or extended for a broader population.

One limitation of this study is that the impact of age and gender on T–E values was not investigated. Table 1 shows that the two groups (CAN− and CAN+) are not completely age matched (mean age difference is 9 years); however, the age difference was found to be not significant (p > 0.05) using the nonparametric Kruskal–Wallis test. In the CAN+ group gender it is equally distributed. However, in the CAN− group, the number of male subjects is higher than the female participants. Since, it is not a completely skewed distribution of gender in both groups and the T–E analysis has achieved 100 % classification result, we may deduce that the result is minimally affected by the gender distribution.

6 Conclusion

Alteration of autonomic nervous system function measured by multi-lag T–E can be used to identify CAN subjects with higher accuracy than conventional T–E method with shorter ECG recording length. Further research on a larger sample size is required to further elucidate the findings of this study and effectiveness of multi-lag T–E analysis for differentiation between mild and definite CAN+ in diabetic patients.