1. Introduction
Sudden Cardiac Death (SCD) has been defined as the unexpected natural and cardiac originated death within a short period of time, given at less than 1 h from the onset of symptoms if witnessed or within 24 h of having been observed alive if unwitnessed, in a person with no prior condition that could be considered as fatal [
1,
2,
3]. There are other causes of sudden death, but cardiac is the most usual origin, specifically including: (a) Arrhythmic causes, such as ventricular tachycardia (VT), ventricular fibrillation (VF), or asystole; (b) several other structural heart disease origins, for instance, those ones corresponding to congenital heart disease; and (c) abnormal function of the autonomous nervous system (ANS), which is not itself a death cause, but it can promote others such as arrhythmic or hypertensive death [
4]. The SCD mechanism in the last case is often VT or VF. Given the incidence of SCD as major cause of mortality in the world, methods have been proposed aiming to provide risk stratification tools for cardiac patients [
5]. SCD episodes can happen not only in patients with coronary or cardiomyopathic disease, but also they can occur in people with no previous heart alteration, which makes the risk stratification extremely complex. The prognostic significance of noninvasive studies and the efficacy of the therapeutic actions have been pointed to be etiology dependent [
6]. The most widely used SCD-risk marker in the clinical practice is the Left Ventricular Ejection Fraction (LVEF), but given its low specificity, many other techniques have been proposed. A relevant subset of them is given by the computational indices that are obtained from the signal analysis of the Electrocardiogram (ECG), including a variety of proposed biomarkers such as late potentials, heart rate variability (HRV), T-wave alternans, or deceleration capacity. The interested reader can refer to [
7] for a detailed review on issues related with signal processing, technology transfer, and scientific evidence for all of them.
Many of these SCD markers are obtained in a Holter recording, which is a diagnostic tool consisting of 24 to 48 h signal registers in two or three chest leads to be subsequently processed by using a computer program, so that a variety of cardiac events can be identified by the clinician. Probably one of the most scrutinized markers of SCD risk from Holter recordings is HRV, which measures the time changes between consecutive cardiac beats [
8]. Its interest partially comes from its non-invasive nature and its easy for analysis, only needing to know the time instants of the beat occurrences. The heart does not behave like a periodic oscillator, but instead its rhythm is modulated by the ANS, and the simultaneous actuation of its two branches (sympathetic and parasympathetic) causes dynamic oscillations of the cardiac frequency, producing the presence of HRV [
9]. Among the many methods that have been proposed in the literature to quantify the HRV indices, the nonlinear methods extract relevant information from HRV signals in terms of their complexity. Nonlinear indices are based on the underlying idea that fluctuations in the between-beat intervals (also known as
intervals) can exhibit characteristics that are well known from Complex Dynamic Systems Theory, and broadly speaking, healthy states are expected to correspond to more complex patterns than pathological states. However, some pathologies are associated with highly erratic fluctuations with statistical properties resembling uncorrelated noise [
7], and traditional algorithms could yield higher irregularity indices for such pathological signals when compared to healthy dynamics, even though the latter represent more physiologically complex states [
10]. This possible inconsistency may be due to the fact that traditional algorithms are based on single scale analysis, and they can not take into account the complex temporal fluctuations inherent to healthy physiologic control systems. It is usual that studies based on 1-day Holter monitoring [
11] envision that relevant information could be obtained from longer duration recordings, however, few studies [
12,
13] have scrutinized nonlinear indices in several-day Holter monitoring, despite its current and increasing availability in the clinical practice. Note in the following that, whereas some authors refer to long-term Holter as those with duration about 24 h, we will use long-term to refer to the Holter recordings when measured for several days throughout this paper.
Several of the nonlinear HRV measurements (based either on Chaos Theory, Information Theory, or Fractal Theory) have been paid special attention according to the electrophysiological hypothesis that the long-term regulation is a homeostatic yet dynamical equilibrium, which can be expected to be complex and multi-cause enough to require a set of indices that should be calculated at different scales. This has motivated the extension of several of those indices to what can be called their
multiscale versions. Remarkable examples of this effort are the Multiscale Entropy (MSE) method [
14,
15,
16,
17], the Multiscale Time Irreversibility (MTI) method [
15], and the Multifractal Spectrum (MFS) method [
18,
19]. MSE has been applied to predict stroke-in-evolution in acute ischemic stroke patients using one-hour ECG signals during 24 h [
20], and also recently to predict vagus-nerve stimulation outcome in patients with drug-resistance epilepsy, who were found to have lower preoperative HRV than controls [
21]. Different preprocessing stages have been proposed for it, including the use of MSE in the first difference of RR-interval time series instead of the series itself, yielding better statistical support and discrimination capabilities between Congestive Heart Failure (CHF) and control groups [
22]. Some authors consider the MSE algorithm biased for two reasons: First, the similarity criteria is fixed for all scales, whereas coarse grained time series variance has been pointed to decrease with the scale; And second, spurious oscillations are introduced due to the suboptimal procedure for eliminating the fast temporal scales of the time series. Accordingly, a modified algorithm was proposed [
23], so-called the refined MSE (RMSE), in order to overcome these limitations, and it was tested in simulations and in 24-h HRV signals from aortic stenosis and control groups. The use of RMSE did not allow to make inferences that could not be made by MSE with real data, however, simulations showed that RMSE can be a more reliable method for the assessment of entropy-based irregularity. A comparative study between MSE and RMSE was performed in [
24] confirming that despite the differences they both present similar tendencies with scale factor. MTI has been applied to HRV and blood preasure variability (BPV) signals concluding that TI of beat-to-beat HRV and BPV is significantly altered during orthostasis [
25]. In addition, recent interest has been raised in the use of some of the multiscale indices in the analysis of atrial fibrillation (AF) dynamics, which has been scrutinized in the context of ischemic stroke prediction in patients with permanent AF [
26].
We can say that cardiac long-term monitoring (LTM) has been technologically achieved, and previous studies exist which have scrutinized the value of these recordings in simple and well-known clinical indices from practice, such as the number or rate of ectopic beats or the number or rate of different cardiac events [
27]. However, algorithm robustness should be paid attention if deeper physiological and pathological information is to be extracted from nonlinear multiscale indices, in order to be sure about their reliability when working with long series in populational data, and to our best knowledge, few previous works can be found noting this point with attention. Therefore, we propose here to study the robustness of nonlinear multiscale HRV measurements to characterize different cardiac health states in LTM recordings.
Specifically, we made a comparative analysis of the three mentioned multiscale methods on a database consisting of patients with 7-day Holter in two cardiac conditions, namely, CHF and AF [
28]. This study aims to give basic knowledge on the usefulness and current limitations of these methods towards their future and principled use for SCD risk stratification. For this purpose, a nonparametric statistical test is proposed in order to compare and give cut-off comparison levels between two different situations in terms of the confidence intervals (CI) for the median difference across multiscale representations of poblational representations. This method can be used either for establishing comparisons among patients or subjects with different conditions, or to scrutinize the impact of preprocessing or data length on the statistical properties of the multiscale representations. The procedure can be seen as an extension of previously used statistical comparisons [
20,
21] in terms of nonparametric bootstrap tests for confidence bands [
12,
29].
The structure of the paper is as follows. In the next section, the fundamentals of the multiscale methods selected here for HRV analysis are described. Then, existing methods and the new procedure based on nonparametric bootstrap tests for median difference are provided in
Section 3, together with the presentation of the available recordings in Physionet [
11] used as starting benchmark for this study and the LTM-ECG databases for CHF and AF during 7 days. In
Section 4, a set of experiments are conducted and results are presented on the suitability, together with some technical limitations and consistency properties, of these benchmarked multiscale algorithms. Finally, in
Section 5, technical and clinical discussion is given and conclusions are drawn.
2. Fundamentals of Multiscale Methods for HRV Analysis
HRV measurements aim to give a numerical magnitude of the time fluctuations between sets of consecutive beats. The short-term recordings of HRV are usually measured about 3 to 5 min, and they have been traditionally associated with the dynamic control of the ANS on the heart rate and the cardiac properties. The long-term fluctuations of HRV have been described to have a wide physiological meaning in terms of the cardiovascular system self-regulation mechanism description. The ANS is divided into two branches, namely, the sympathetic and the vagal (parasympathetic) ones. Broadly speaking, the activation and excitation of the sympathetic branch has an accelerating effect on the cardiac cycle, whereas the vagal activation has a decelerating effect, but both subsystems are simultaneously and continuously working and compensating themselves, so that oscillations on a dynamic equilibrium are produced on the heart rate [
9]. In addition, the ANS receives information through the so-called efferent pathways from a wide variety of systems and organs (heart, digestive system, kidney, respiratory system, and many others), and those influences are part of the genesis of ANS afferent pathways, in which heart rate is affected and is involved through different control mechanisms. Additional influences on the ANS such as humoral factors, night-day cycles, or environmental influences, are slower than the ANS, so that they only can influence the long-term HRV [
30]. All these sources are contributors to the HRV modulation, which globally has its origin on a complex dynamic equilibrium arising from diverse mechanisms in the cardiovascular system that are taking place in the short-, the middle-, and the long-term scales.
A number of scientific, technical, and medical studies have focused on the HRV, and it well might be the most studied index in the SCD risk-stratification literature. Nonlinear methods [
9] include several subfamilies, according to their calculation being based on Information Theory, on Chaos Theory, or on Fractal Theory. These methods have been paid special attention not only for their attractive theoretical foundations, but also because they seemed to have promising risk capabilities in small-sized studies with patients [
7]. A usual situation in the literature of nonlinear HRV indices has been that some basic index has been first proposed, which has shown descriptive capabilities and some independence for SCD risk stratification, and then this index has been subsequently extended to a multiscale formulation, aiming to capture a richer variety of descriptions for the signal behavior. As described next, this is the case of MSE, MTI, and MFS methods for HRV analysis.
We denote the continuous-time ECG signal of a patient by
, and we register it during an observed time interval, denoted by
. As a result of preprocessing steps devoted to signal filtering and R-wave detection, we can detect the R wave of each beat in
, and the time instants associated to each R-wave are denoted as
, with
, so that the detected set of R-waves can be expressed as a point process, given by
where
denotes the continuous-time Dirac’s delta function. It is often useful to work with the so-called
normal beats, which correspond to R-waves in beats that have been only originated in sinus rhythm conditions and where artifacts and ectopic beats are discarded. Here we assume that the R-waves correspond to cardiac beats but not to artifacts or wrong R-wave detections. The RR-tachogram [
9], or just
tachogram, can be denoted by
and is defined as the discrete-time series given by the indexed time difference between consecutive R-wave times (excluding artifacts and non-physiological beats, as well as conventional quality-control beat filters), this is,
The following algorithms and indices can be expressed and obtained in terms of the tachogram registered in patients.
2.1. MSE Analysis
The approximate entropy (
) can be described as a nonlinear fluctuation measurement that aims to quantify the irregularity of a RR-interval time series [
31]. An
increase is usually interpreted as an indicator of irregularity increase in the underlying cardiovascular process. The Sample Entropy (
) was subsequently introduced [
32] to solve the limitations of the
[
33], since the latter compares each pattern in a time series with other patterns but also with itself, which leads to the overestimation of similarity existence in that time series, and hence to strong bias and to some inconsistent results.
is the negative of the natural logarithm of the conditional probability that two similar patterns of
m point segments,
and
of tachogram
, remain similar if we increase the number of points to
, within a tolerance
r that is defined as a noise-rejection filter [
32]. SampEn index reduces the statistical bias of ApEn index and it is a measurement rather independent of the data length, but its problem is that it it can be unstable when the counted events are scattered.
The new concept of multiscale analysis was proposed to overcome several limitations pointed out for
and
measurements which could be leading to clinical misinterpretation of HRV in some conditions. The MSE analysis was introduced by Costa et al. in [
14], and it was also intended to provide a richer description of the cardiac dynamics in terms of a set of naturally related indices, rather than to use a single number. For a given discrete-time series, a new series is constructed in a scale
, the terms of which are the average of the consecutive elements of the original series without overlapping. For a time series with
, this corresponds to the original series, whereas for
, the series is constructed with the average of the elements taken from two by two, and so on. We finally calculate the
for each one of these new generated series. When the obtained values are represented versus the scale factor, the dependence of the measured entropy with the time scale can be scrutinized. The maximum scale to use depends on the number of samples in the time series.
Starting with tachogram signal , we denote as to explicitly consider the dependence of its design parameters. We obtain the consecutive time series , determined by scaling factor as follows:
First, the original time series is divided into non-overlapping intervals with window size of samples. Then the signal mean is obtained for each of the sample windows.
Each element of the series
is calculated according to the equation:
where
for notation simplicity. For the first scale, the time series
is just the original time series. The length of each obtained time subseries is equal to
.
The sample entropy index is calculated for each time series , and it is represented as a scale-factor function .
The MSE analysis has been applied to a variety of cardiac and cardiopathological situations, including the analysis of CHF [
14], hypertensive and sino-aortic denervated conditions in experimental studies [
24], and the ANS evolution before, during, and after percutaneous transluminal coronary angioplasty [
34], as well as for prediction of ischemic stroke in patients with persistent AF [
26], among other things. In all these cases, the length of the analyzed signals did not exceed 24 h using 20 as a maximum scale value, which involves
tachograms with about or more than 20,000 samples each.
2.2. MTI Analysis
Time irreversibility has recently attracted attention in the cardiovascular-signal field. A signal is said to be time irreversible if its statistical properties change after its time reversal. The consistency loss of the statistical properties of a signal when the signal reading undergoes a change through time inversion is measured using the MTI, which represents an asymmetry index. This index is higher in healthy systems (with more complex dynamics) and it decreases in conditions like pathology or aging, as introduced in [
15,
16]. On the other hand, physiological time series generate complex fluctuations in multiple depending time scales, due to the existence of different hierarchical and interrelated regulatory systems.
The MTI analysis has been applied to measure nonlinear dynamics in heart-rate time series. For instance, MTI indices were computed in [
35] for 20 healthy neonates to detect the presence of nonlinearity in their cardiac-rhythm control system, and temporal asymmetries were detected within their heart rate dynamics even shortly after birth.
2.3. MFS Analysis
Physiological signals have been shown to present fractal temporal structure under healthy normal conditions [
12,
36,
37]. In particular, it was shown in [
38] that time series generated by certain cardiovascular control systems in healthy conditions require a large number of exponents to adequately characterize their scaling properties, and that the nonlinear properties of this behavior are encoded in the Fourier phases. The same work used examples of CHF patients to contrast the previous finding with the loss of multifractality in this example of life-threatening condition. In this setting, RR-interval time series have been analyzed in terms of multiple scaling exponents. We next summarize the basic principles for estimating the MFS from RR signals that are usually followed in the literature. The interested reader can consult the original works on its application [
19,
38] and the details on the wavelet-transform modulus maxima method [
39], which gives a principled calculation method for this purpose.
Whereas monofractal signals have the same scaling properties through time and they can be indexed by a single global exponent (such as Hurst exponent
H) characterizing their fluctuations, other signals exhibit variations in their local Hurst exponent along time. When several subsets of a signal are characterized by the same local Hurst exponent
, and when each of these signals can be characterized by a fractal dimension measurement, we denote this estimated dimension as
. Accordingly,
will have nonzero values on a set of discrete points in
h for some class of signals. Local value of
h is modernly estimated with Wavelet Theory [
39], often using successive derivatives of the Gaussian function as the analyzing wavelet at different scales
a, in order to remove polinomial trends with polinomial order up to the wavelet derivative order. In these conditions, the problem reduces to obtain the modulus of the maxima extrema of the time series wavelt transform at eath time instant.
We then estimate the partition function
, as the sumation of the
qth powers of these local maxima as a function of scale, and for small scales, it is fulfilled that it has the form
, where
are exponents that can be estimated. In monofractal signals, a linear scaling exponent spectrum is obtained, given by
. However, for multifractal signals we obtain a nonlinear expression, and it can be shown that
, where
is not constant. Accordingly, the estimated fractal dimensions
are obtained by the Legendre transform of
, finally yielding
Given that can be related to uncorrelated changing time series, this representation allows to determine to what extent a process conveys anticorrelated () or correlated () behaviour consistently manifested through different time periods.
The multifractal structure of HRV can reflect important properties of the heart-rate autonomic regulation. Multifractal analysis is an expansion of fractal analysis since it characterizes the time series variability with a collection of scaling exponents instead of a single one, which makes possible to investigate and quantify HRV in terms of its multiexponent properties. A right shift has been revealed in the multifractal spectrum peaks for healthy subjects during meditation [
40], which points to a better health condition of persons with respect to multifractal nature. Accordingly, a healthy heart-rate regulation promotes a multifractal signal.
5. Discussion and Conclusions
In this work, we have addressed the calculation of three representative multiscale indices, namely, MSE, MTI, and MFS, on 1-day and 7-day Holter recordings. From our results, we can conclude that, when present, the trends are consistent between the additional scales provided by 1-day recordings and by the 7-day recordings, but the second ones can give a statistically better view of this kind of representations, specially in the MFS representations.
Contributions. Several preceding results can be found in which descriptions of 1-day estimated indices are shown, however, and to the best knowledge of the authors, few works address from this perspective the 7-days case with scales up to 100. Accordingly, and given the advances in the monitoring and well-being technology in our days, this represents a good moment to determine the possible advantages and issues of using these methods in larger observation scales. Non-linear methods used in this work (MSE, MTI, and MFS) have previously been proposed and widely used in the literature. The contributions of our work are twofold. On the one hand, we analyzed the dynamics, the consistency, and especially the robustness of these indices when applied to the new long term scenarios (7 days) and to point out that algorithms may need to be improved when using them in this area. On the other hand, to find out if further information is conveyed in the temporal scales that have not been analyzed before in shorter recordings.
Summary and Discussion of Results. Results showed significant differences between CHF and AF populations not only for short-term scales but also for long-term and very long-term scales in some indices, MSE (MTI) being higher (lower) for AF. These results are consistent with previous studies [
14,
15,
16]. Here we confirmed the same trends, but we also obtained several differences in very long-term scales that had not been analyzed before. MSE was not further significantly different for very long-term scales. This may be attributed to the fact that statistical characteristics of AF signals resemble those from noise, and being MSE a entropy-based measure, it is higher for AF than for CHF for short-term and long-term scales, but not for very long-term scales, where the surrogate signals present attenuated these erratic characteristics. On the other hand, MTI shows that time irreversiblity is for all time scales lower in AF patients, i.e., for those presenting a more severe pathology.
We detected that 1-day estimations of the multiscale indices can show a distorted profile. The different estimations observed in the different time windows could be due to several reasons. In addition, this effect could have statistical roots, in terms of consistency and variance of the estimation with shorter time series. On the other hand, it also could be due to the fact that some frames had some properties in their dynamics which made the multiscale algorithm fail or disrupt. The differences in daily activities should not be a limitation, as far as the multiscale algorithm measurements are the correlations present in the signal at different time scales. In the conventionally used surrogate-signal test, which is built with lower-pass versions of the original signal, the more we increased the scales, the more the short-term relationships were eliminated to manifest the long-term. Thus, the short-term differences would be hidden for high scales. We thoroughly checked that the underlying signals did not have artifacts or aberrations with a specifically customized software to support an observer to watch at the heart-rate signals and at the original ECG signals simultaneously. Accordingly, one of our conclusions is that these multiscale algorithms, while likely informative, should trust not only on the increase of the length of the signals for their consistency, but also on improving their algorithmic robustness. While stationarity can not often be assumed for HRV signals, MSE was originally proposed and applied in 24-h Holter recordings, and later many other works have used MSE to study different dynamics in 24-h. The analysis of the underlying dynamics during the day or night can be complex with these methods. Nevertheless, we observed that the disruption effect in one day does not have any significant impact in the 7-day analysis, which promotes their use in favor of improved statistical consistency.
Usefulness of Median-difference Tests. An extension of previously proposed statistical comparisons for scales has been proposed here in terms of nonparametric bootstrap resampling, which allows us to establish poblational comparisons in terms of the median difference of the multiscale representations, either for different health conditions, or for different acquisition conditions. Our results with this method were consistent with the results identified in the literature and showed some non-observed differences, specially the ones related with AF patients, and with different statistical consistency behavior and retrieved information about the patient given by these three multiscale indices, which suggests their use jointly in these and other populations.
Related Relevant LTM Applications With the technological availability of cardiac monitoring systems for extended time periods, LTM is expected to bring interest to new applications in health, wellness, and research. For instance, it has been recently pointed out [
30] that energetic environmental phenomena can affect psychophysical processes on people in different ways depending on their sensitivity, health status, and ability of the ANS to self-regulate. In that study, the HRV was recorded for 72 consecutive hours per week over a five-month period in 16 participants, in order to examine ANS responses during normal background environmental periods. Interestingly, HRV measurements were negatively correlated with the solar wind speed, and the low-frequency and high-frequency power were negatively correlated with the magnetic field. This study confirmed that the daily activity of the ANS responds to changes in geomagnetic and solar activity during periods of undisturbed normal activity, it starts at different times after changes in various environmental factors, and it persist for variable time periods. As far as the activity of the ANS reflected by HRV measurements is affected by solar and geomagnetic influences, the analysis of HRV should take into account these effects when possible. This study is focused on the ANS modulation, and hence spectral measurements were mostly used, nevertheless, this kind of data could be analyzed with multiscale indices to provide a wider view on the long-term behaviour of HRV.
On the Clinical Usefulness of Multiscale Indices. As described with detail in [
7], many indices have been proposed in the last years to develop risk stratification from different kinds of analysis of electric cardiac signals. However, these techniques from the academic research world rarely are used in the clinical practice. In that reference, our team analyzes the possible reasons by decoupling the sometimes limited accuracy and the lack of consensus on the robustness with the appropriate signal processing implementations. In this line of search, the present work aims to first establish the need for robustness in the methods that are widely used nowadays, as a requirement before enrolling in risk stratification studies, which require high-cost and high effort to yield clinically useful use to these techniques. Our future research would consider first to perform the same multiscale analysis (MSE, MTI, and MFS) in LTM healthy subjects, in order to study the dynamic behavior in normal conditions when increasing the number of scales. As indicated, the development of more robust multiscale indices is a desirable target in order to continue to progress in this informative characterization of the cardiac dynamics. On the other hand, AF and CHF are different heart diseases with different dynamics, though in this study we only could analyze AF patients as a subset of CHF patients. CHF is a syndrome of the deterioration in short-term and long-term regulatory mechanisms, while AF, and especially the persistent one, is an intrinsic short-term mechanism. Whereas long-term mechanisms can be present in AF, it seems that LTM recordings should be better used to analyze their presence or deterioration.
Despite the need for more robust algorithms in long-term nonlinear indices, we still consider that it worth the effort to use these indices, as vast literature supports their informativeness in other scenarios in addition to the very long-term monitoring. In general terms, the underlying hypothesis in preceding studies is that heart rate oscillations respond to phenomena with very different characteristic times, ranging from seconds to months. The former are well evaluated in a 24-h recording, but the low frequency (such as circadian cycles, many hormonal cycles, or secondary to changes in activity during the week) would only be represented in long-term monitoring. The response to these stimuli has been moderately studied, but it could provide us with interesting information on some clinical aspects, such as the prediction of decompensation in heart failure, the evolution of cardiovascular remodelling after acute injury (e.g., in the weeks following acute myocardial infarction), the evolution of sports training, or the susceptibility to the development of malignant arrhythmias (related with sudden cardiac death). All these aspects involve hormonal processes, inflammatory processes, adaptations of different organs or systems, and they are characterized by the interaction of mechanisms with slow response times, so that they develop in days or weeks, rather than in minutes or hours. Therefore, methods of analyzing the very long-term behavior of heart rate responses could be of interest, especially in noisy rhythms such as AF, in which almost all the usual HRV parameters are difficult to interpret.
We still do not know if these methods will be widely used in future clinical practice because we are still in the phase of describing how these indices behave, their variability and their comparison with the results in 24 h, among other thing. In any case, all these efforts should be supported by robust algorithms for multiscale nonlinear indices, an idea that had not been previously paid much attention in the nonlinear HRV literature.