K-nearest-neighbor conditional entropy approach for the assessment of the short-term complexity of cardiovascular control

Given the stationary time series y = {y(i), i = 1, ..., N} let us define as a pattern of length L the ordered sequence of L samples, y_L(i) = (y(i), y(i − 1),..., y(i − L + 1)). The pattern y_L(i) is actually a point in the L-dimensional embedding space reconstructed with the technique of the lagged coordinates with delay equal to 1. The pattern y_L(i) can be seen as the sequence formed by the current sample, y(i), and by the sequence of L − 1 past samples, y_L−1(i − 1). The sample y(i) is referred to as the image of y_L−1(i − 1) in the following. We define as y_L = {y_L(i), i = L, ..., N} and y_{L − 1} = {y_{L − 1}(i − 1), i = L, ..., N} the sets of patterns of length L and L − 1 respectively. The CE associated with y_L measures the average amount of information carried by the most recent sample of a pattern in y_L, y(i), when L − 1 previous samples (i.e. y_L−1(i − 1)) are given. It is computed as the weighted sum of the Shannon entropy (SE) of the images of y_L−1(i − 1), y/y_L−1(i − 1):

$$\begin{equation} {\rm CE}(L) = \sum {p(y_{L - 1} (i - 1)) \cdot {\rm SE}(y/y_{L - 1} (i - 1))} , \end{equation} \tag{ 1 }$$

where the weight is the probability of the pattern y_L−1(i − 1), p(y_L−1(i − 1)). The probability p(y_L−1(i − 1)) is estimated as the fraction of patterns in y_L−1 similar to y_L−1(i − 1). The SE of the images of y_L−1(i − 1) is given by

$$\begin{equation} \fl {\rm SE}(y/y_{L - 1} (i - 1)) = - \sum {p(y(i)/y_{L - 1} (i - 1)) \cdot \log (p(y(i)/y_{L - 1} (i - 1))} ), \end{equation} \tag{ 2 }$$

where log is the natural logarithm and

$$\begin{equation} p(y(i)/y_{L - 1} (i - 1)) = p(y_L (i))/p(y_{L - 1} (i - 1)) \end{equation} \tag{ 3 }$$

is the conditional probability that the current sample assumes the value y(i) given that its L − 1 previous samples are y_L−1(i − 1). While the sum in (1) is extended over all the possible patterns y_L−1(i − 1), the sum in (2) is extended over the possible values y(i) given y_L−1(i − 1). CE is bounded between 0 and SE of y. CE = 0 is computed when the distribution of the images of y_L−1(i − 1) exhibits a unit peak and this situation occurs for any y_L−1(i − 1), thus suggesting the full predictability of the current sample given L − 1 past values. The SE of y

$$\begin{equation} {\rm SE}(y) = - \sum {p(y(i)) \cdot \log (p(y(i)} )) \end{equation} \tag{ 4 }$$

represents the average information carried by y(i) when no previous samples are given. In addition, CE(L) is always decreasing (or constant) with L since the increase of the number of past values reduces (or leaves unchanged) the uncertainty about future values. CE(1) is set to SE(y).

The CE is estimated according to two strategies.

The first strategy relies on the construction of the distribution of the images of y_L−1(i − 1), y/y_L−1(i − 1), on the estimate of SE(y/y_L−1(i − 1)) and on the computation of (1). The construction of y/y_L−1(i − 1) implies scanning y_L−1 to search for y_L−1(i − 1). When a pattern y_L−1(j − 1) similar to y_L−1(i − 1) is found in y_L−1, the image of y_L−1(j − 1), y(j), is stored, thus constructing y/y_L−1(i − 1). Given y/y_L−1(i − 1), SE(y/y_L−1(i − 1)) can be computed.

The second strategy does not need the estimation of SE(y/y_L−1(i − 1)). It is based on the direct estimation of the conditional probability of y(i) given y_L−1(i − 1), p(y(i)/y_L−1(i − 1)), according to (3). It implies scanning y_L and y_L−1 to search for patterns similar to y_L(i) and y_L−1(i − 1) respectively. The ratio of the fraction of patterns similar to y_L(i) in y_L to that of patterns similar to y_L−1(i − 1) in y_L−1 allows the estimation of p(y(i)/y_L−1(i − 1)).

Three traditional approaches to the estimate of entropy-based complexity will be described in the following. One of these approaches makes use of the first strategy (i.e. the UQ approach), while the remaining two methods exploit the second strategy (i.e. ApEn and SampEn).

This approach makes use of the first strategy described in section 2.2. The full dynamics of the series is spread over ξ bins of size ε = (max – min)/ξ where max and min stand for the maximum and the minimum of y. The UQ imposes a perfect partition of the L-dimensional embedding space into ξ^L cells (i.e. hyper-cubes of side ε). Cells are disjoint and their union covers the entire L-dimensional embedding space. Cells define the coarse graining of the embedding space and similarity criterion: indeed, all the patterns inside a cell are indistinguishable. Therefore, given a pattern y_L−1(i − 1), every pattern belonging to the same cell of y_L−1(i − 1) contributes with its image to form the distribution of y/y_L−1(i − 1). Given the UQ procedure the distribution of y/y_L−1(i − 1) is quantized over ξ bins and SE(y/y_L–1(i − 1)) is easily computed from the sample frequencies. The evaluation of SE(y/y_L−1(i − 1)), extended to all the patterns y_L−1(i − 1), allows the assessment of the CE according to (1).

Approximate entropy (ApEn) makes use of the second strategy described in section 2.2. Indeed, ApEn is computed as

$$\begin{equation} {\rm ApEn}(L,r) = - \frac{1}{{N - L + 1}}\sum\limits_{i = L}^N {\log \frac{{N_i (L,r)}}{{N_i (L - 1,r)}}} , \end{equation} \tag{ 5 }$$

where N_i(L, r) is the number of points in y_L at distance smaller than r from y_L(i) and N_i(L − 1, r) is the number of points in y_L−1 at distance smaller than r from y_L−1(i − 1). As a consequence ApEn is calculated without approximating SE(y/y_L−1(i − 1)) and solely based on the estimate of (3) according to the fraction of patterns similar to y_L−1(i − 1) in y_L−1 that remain similar, within a tolerance r, in y_L. Since we utilize the Euclidean norm to measure the distance between patterns in the embedding space, cells defining coarse graining of the embedding space are hyper-spheres of radius r, constructed around a given pattern (i.e. N − L + 1 intersecting pattern-centered cells). At difference with the UQ approach, cells intersect each other. However, all the cells have the same dimension as in the UQ approach. It is worth stressing that since y_L(i) is always at distance closer than r from itself, both N_i(L, r) ≥ 1 and N_i(L − 1, r) ≥ 1. In addition, since N_i(L − 1, r) ≥ N_i(L, r), N_i(L, r)/N_i(L − 1, r) ≤ 1. ApEn(1) is computed as the negative average natural logarithm of the probability that two samples are at distance closer than r.

Similarly to ApEn, SampEn is based on the direct approximation of the conditional probability. However, instead of approximating (3) as in ApEn, SampEn estimates the conditional probability that two patterns that are similar in y_L−1 remain similar, within a tolerance r, in y_L. It is given by the ratio of ∑^N_{i = L}N_i(L, r) to ∑^N_{i = L}N_i(L − 1, r). Therefore, SampEn is defined as

$$\begin{equation} {\rm SampEn}(L,r) = - \log \frac{{\sum\nolimits_{i = L}^N {N_i (L,r)} }}{{\sum\nolimits_{i = L}^N {N_i (L - 1,r)} }}. \end{equation} \tag{ 6 }$$

The Euclidean norm is utilized as in ApEn and, consequently, the coarse graining procedure is the same (i.e. N − L + 1 intersecting pattern-centered hyper-spheres). Since N_i(L − 1, r) ≥ N_i(L, r), then ∑^N_{i = L}N_i(L, r)/∑^N_{i = L}N_i(L − 1, r) ≤ 1. SampEn(1) is computed as the negative natural logarithm of the probability that two samples are at distance closer than r.

The UQ approach, ApEn and SampEn share a common feature: the dimension of the cells for coarse graining the embedding space is fixed (Porta et al 2007c). As a consequence, when the time series are relatively short, these techniques cannot ensure a reliable statistic for the estimation of the conditional distribution because very few patterns might be present inside the cells. Moreover, the reliability of the estimate deteriorates more and more with the length of the conditioning pattern (i.e. L − 1), thus leading to completely unreliable conditional distributions at high L. If the series are short, the complete loss of the reliability occurs at very low L (Porta et al 1998). The full unreliability occurs as early as only one pattern is present in the cell. This situation is commonly referred to as 'self-matching' in ApEn and SampEn calculation (Richman and Moorman 2000). In ApEn and SampEn calculation this unique pattern is located in the center of the cell, while in the UQ approach it can be anywhere inside the cell. 'Self-matching' introduces a significant bias in the estimate of conditional probabilities and, consequently, on the estimate of CE. Indeed, p(y(i)/y_L−1(i − 1)) = 1 given the unique appearance of y_L−1(i − 1). As a result the logarithm is equal to 0, thus artificially reducing CE. Therefore, strategies have been implemented to reduce this bias. The rationale of these strategies is mainly to substitute the false certainly associated with the unique appearance of the pattern with the maximum uncertainty that can be estimated over the series. In the case of the UQ technique the CE has been corrected, thus defining the CCE, by adding to a term proportional to the fraction of 'self-matching' patterns (Porta et al 1998). The coefficient of proportionality is the SE(y), thus substituting the erroneous null contribution of the 'self-matching' pattern to the CE estimate with the maximal amount of information carried by y (Porta et al 1998). In the case of ApEn a corrected ApEn (CApEn) has been defined by substituting the ratio N_i(L, r)/N_i(L − 1, r) with 1/(N − L + 1) when N_i(L, r) = 1 or N_i(L − 1, r) = 1 (Porta et al 2007b), thus forcing a situation of maximum uncertainty in place of false certainty. In the case of SampEn, correction is simply obtained by excluding 'self-matching' patterns from the computation of N_i(L, r) and N_i(L − 1, r) (i.e. both N_i(L, r) and N_i(L − 1, r) are diminished by 1) (Richman and Moorman 2000).

In order to overcome the issue of 'self-matching', a completely different strategy can be followed. Instead of exploiting a coarse graining procedure based on equal size cells, we make use of a coarse graining approach based on cells of different sizes (Porta et al 2007c). This approach makes use of the Euclidean norm to measure the distance between patterns and the same coarse graining procedure as ApEn and SampEn (i.e. N − L + 1 intersecting pattern-centered cells). However, instead of keeping constant the size of the cells, it is adjusted to include a fixed amount of points (i.e. the k nearest neighbors of y_L−1(i − 1)). All the k nearest neighbors of y_L−1(i − 1) contribute with their images to form the conditional distribution y/y_L−1(i − 1), thus limiting the loss of reliability of y/y_L−1(i − 1) with L. SE(y/y_L−1(i − 1)) is easily computed after UQ of y or as the negative natural logarithm of the probability that two samples of y/y_L−1(i − 1) are at distance closer than r as in SampEn. CE is estimated as

$$\begin{equation} {\rm KNNCE}(L) = \frac{1}{{N - L + 1}}\sum\limits_{i = L}^N {{\rm SE}(y/y_{L - 1} (i - 1))} . \end{equation} \tag{ 7 }$$

In addition, this approach shares with CCE a significant advantage: regardless of the dynamics of the process, k-nearest-neighbor CE (KNNCE) exhibits a well-defined minimum over L. Since the appendage of one dimension to the embedding space (i.e. the unit increase of the length of the pattern) causes the spread of all the points in the embedding space, the k nearest neighbors of y_L−1(i − 1) move apart. The spread produces two distinct consequences: (i) the unfolding of the dynamics, thus resolving the ambiguities in fixing future behaviors and decreasing KNNCE; (ii) the increase of uncertainty about future values, thus weakening prediction and increasing KNNCE. The balance between these two opposite contributions leads to the formation of a KNNCE minimum over L. Even though a minimum is always detectable even in the presence of fully unpredictable dynamics (e.g. white noise), its deepness depends on the type of dynamics, being deeper when repetitive patterns are found.

This protocol was originally designed to study the effects of the pharmacological blockades of the parasympathetic and sympathetic branches of the autonomic nervous system on baroreflex sensitivity (Parlow et al 1995, Toader et al 2008). Briefly, we studied nine healthy male physicians aged from 25 to 46 years familiar with the study setting. None of them had any abnormal finding in history, physical examination or electrocardiography or was receiving any medication. All had normal resting brachial arterial pressure measured with a sphygmomanometer. They were instructed to avoid tobacco, alcohol and caffeine for 12 h and strenuous exercise for 24 h before each experiment. Electrocardiogram and noninvasive finger blood pressure (Finapress 2300, Ohmeda, Englewood, Colorado, USA) were recorded during the experiments. The hand of the subject was kept at the level of the heart. Signals were sampled at 500 Hz. Experimental sessions were performed in 3 days at approximately 2 week intervals. One volunteer took part only in the first day experiments. Subjects remained at rest in supine position in a quiet darkened room during all the recordings. Each experiment started in the morning between 08:00 and 09:00 AM and consisted of 15–20 min of baseline (B) recording followed by 15–20 min of recording after drug administration. Recordings were obtained: (i) on day 1 after parasympathetic blockade with 40 µg kg⁻¹ i.v. atropine sulfate (atropine, AT) to block muscarinic receptors; (ii) on day 2 after β-adrenergic blockade with 200 µg kg⁻¹ i.v. propranolol (PR) to block β₁ cardiac and β₂ vascular peripheral adrenergic receptors; (iii) on day 1 PR was administered at the end of the AT session (the dose of AT was reinforced by 10 µg kg⁻¹) to combine the effect of AT and PR (AT + PR) and obtain a cardiac parasympathetic and sympathetic blockade; (iv) on day 3 recordings were obtained 120 min after 6 µg kg⁻¹ per os clonidine hydrochloride (clonidine, CL) to centrally block the sympathetic outflow to heart and vasculature. All the subjects gave their written informed consent. The protocol adhered to the principles of the Declaration of Helsinki. The human research and ethical review board of the Hospices Civils de Lyon approved the protocol.

After detecting the QRS complex on the electrocardiogram and locating the R-apex using parabolic interpolation, the temporal distance between two consecutive R parabolic apexes was computed and utilized as an approximation of HP. The maximum of arterial pressure inside HP was defined as systolic arterial pressure (SAP). The occurrences of QRS and SAP peaks were carefully checked to avoid erroneous detections or missed beats. If isolated ectopic beats affected HP and SAP values, these measures were linearly interpolated using the closest values unaffected by ectopic beats. HP and SAP were extracted on a beat-to-beat basis. The series were linearly detrended. Sequences of 256 consecutive measures were randomly selected inside each experimental condition. This length of the series was chosen because short-term analysis of cardiovascular variability series is usually performed over about 5 min recordings or 250 samples (Task Force 1996), thus focusing time scales of mechanisms typically involved in short-term cardiovascular regulation. If evident nonstationarities, such as very slow drifting of the mean or sudden changes of the variance, were present despite the linear detrending, the random selection was carried out again.

The following CIs and normalized CIs (NCIs) were assessed: (i) CI_CCE was the minimum of CCE over L after UQ based on ξ = 6 bins (Porta et al 1998, Porta et al 2007a); (ii) NCI_CCE was the ratio of CI_CCE to SE(y) assessed after UQ based on ξ = 6 bins (Porta et al 2000, Porta et al 2007a); (iii) CI_CApEn was assessed with L − 1 = 2 and r = 20% of the standard deviation (Pincus 1995, Porta et al 2007b); (iv) NCI_CApEn was the ratio of CI_CApEn to ApEn(1) (Porta et al 2007b); (v) CI_SampEn was assessed with L − 1 = 2 and r = 20% of the standard deviation (Richman and Moorman 2000); (vi) NCI_SampEn was the ratio of CI_SampEn to SampEn(1) (Porta et al 2007b); (vii) CI_KNNCE was the minimum of KNNCE over L with k = 30 and UQ based on ξ = 6; (viii) NCI_KNNCE was the ratio of CI_KNNCE to SE(y) assessed after UQ based on ξ = 6 bins.

Kruskal–Wallis one way analysis of variance on ranks was applied to check whether CIs (or NCIs) changed after drug administration (Dunn's test for multiple comparisons versus B). After pooling together CIs (or NCIs) regardless of the experimental condition, Friedman repeated measures analysis of variance on ranks was applied to check differences among medians of CIs (or NCIs). Dunn's test was utilized to check the significance of all pair-wise comparisons. After pooling together CIs (or NCIs) regardless of the experimental condition, linear regression analysis was carried out between any pair of CIs (or NCIs). The Pearson product moment correlation coefficient, r, was calculated. If r was significantly different from 0, a significant linear correlation between indexes was detected. The coefficient of variations was assessed to compare the variability of indexes characterized by different medians. The significance of the difference between coefficients of variation was evaluated (Miller 1991). A p < 0.05 was always considered significant.

The most important features of the KNNCE approach are: (i) to limit the loss of reliability of the conditional distribution with the length of the conditioning pattern; (ii) to allow the assessment of CI and NCI according to an optimization procedure.

One of the main disadvantages of CApEn, SampEn and CCE is related to the loss of reliability of conditional distribution as a function of the embedding dimension. This limitation is the consequence of the utilization of a coarse graining strategy of the embedding space exploiting cells of the same size. CApEn, SampEn and CCE deal with this well-known limitation because it leads to a significant reduction of CE (Porta et al 1998, Richman and Moorman 2000, Porta et al 2007b). SampEn prevents the consideration of unreliable conditional distributions (i.e. it excludes 'self-matching' patterns). CApEn corrects unreliable conditional probabilities according to a priori assumptions (i.e. it substitutes erroneous certainty with full uncertainty). CCE adds a corrective term preventing the erroneous decrease of CE estimate (i.e. it imposes an increase of CE estimate according to a measure of unreliability of CE estimate provided by the fraction of 'self-matching' patterns). The KNNCE provides a CE estimate that does not need to introduce any corrective term or a priori information to deal with the loss of reliability of conditional distribution with the pattern length. Indeed, all the conditional distributions have the same consistency assigned a priori according to the number of nearest neighbors (i.e. k). In addition, consistency is independent of the conditioning pattern, y_L−1(i − 1), and of its length: indeed, the number of nearest neighbors is the same for every pattern and constant with L.

Another relevant disadvantage of CApEn and SampEn is related to the arbitrary selection of the length of the conditioning pattern (i.e. L − 1) to derive CI and NCI. The original studies proposed L − 1 = 2 (Richman and Moorman 2000, Pincus et al 1993). This choice is definitely corrected when the dynamics of the series can be largely predicted using very few past samples and when the number of samples helpful to forecast future behaviors is invariant while changing groups and/or experimental conditions. However, a safer procedure should allow the automatic selection of the pattern length most useful to reduce the uncertainty of future behaviors. For example, if a large number of past samples is necessary to solve the ambiguities of the dynamics, as it occurs in the case of figure 1(a), the arbitrary assumption that few samples are sufficient to make a reliable prediction overestimates the complexity of the dynamics (see figures 1(e), (g)). The KNNCE approach allows the estimate of CI_KNNCE and NCI_KNNCE via a minimization procedure over L. It exploits the spread of the dynamics with L leading to the enlargement of the cell containing the k nearest neighbors of y_L−1(i − 1). If the spread is helpful in unfolding dynamics and reducing the uncertainty in predicting future samples, KNNCE decreases. Conversely, if the spread leads to an increase of unpredictability, KNNCE rises. This situation, in association with the bounded nature of KNNCE between 0 and CE(1) = SE(y), produces an unequivocal KNNCE minimum over L. The minimum can be easily extracted via a fully automatic minimization procedure to derive an optimized CI_KNNCE (or NCI_KNNCE).

It is worth noting that CE of y given a pattern formed by L − 1 past samples can be estimated as the difference between the SE of y_L and the SE of y_L−1 (Porta et al 1998) or, equivalently, as the difference between the SE of y and the mutual information (MI) shared by y and the pattern formed by L − 1 past samples (Kraskov et al 2004). However, in this study neither of the two approaches was followed, thus avoiding the bias arising from errors individually made in the estimation of SE and MI (Kozachenko and Leonenko 1987, Kraskov et al 2004) that unlikely cancel each other in the difference.

CIs spanned significantly different ranges of values. Therefore, we recommend to avoid any comparison between studies exploiting different indexes and to draw index-specific conclusions according to the adopted experimental protocol. The index providing the largest coefficient of variation was CI_CApEn, being significantly larger than CI_CCE and CI_KNNCE in the SAP series. The large coefficient of variation of CI_CApEn might suggest a reduced discriminative statistical power of CI_CApEn compared to the remaining ones. Indeed, a larger number of subjects could be necessary to make the type I error probability, p, smaller than a user-defined significance level (e.g. 0.05). Nevertheless, the presumed limited power of CI_CApEn has no effect in this experimental protocol. Normalization of CIs, leading to NCIs, reduced the variability of indexes over the HP series (i.e. all the coefficients of variation of NCIs were smaller than those of CIs), thus prompting for their use in heart rate variability study to improve discriminative statistical power. This finding did not hold when NCIs were assessed over the SAP series. Since the normalization factor accounts for the SE of the series (Porta et al 2007a), we suggest that SE of HP is more responsible for the dispersion of CI indexes than SE of SAP.

Any CI pair was significantly linearly correlated. This result suggests that it is useless to report results derived from all the methods since they do not provide independent information. Preference should be given to the method allowing the clearest distinction between different conditions. However, given that coefficients of variation were not significantly different, except for CApEn, this selection cannot be made a priori. The correlation coefficient was higher between CI_CCE and CI_KNNCE and between CI_CApEn and CI_SampEn. This finding suggests that the largest correlation can be found between methods utilizing the same technique for estimating CE: indeed, CCE and KNNCE computed (1) via the estimation of SE(y/y_L−1(i − 1)), while CApEn and SampEn were based on the straight computation of the conditional probabilities. Normalization generally improved the correlation coefficient, especially when indexes were assessed over the SAP series, thus indicating that the different estimates of SE of the series were responsible for a certain degree of uncorrelation between CIs.

At B in healthy subjects the cardiac pacemaker is under control of vagal and sympathetic branches of the autonomic nervous system (Task Force 1996). The parasympathetic blockade induced by high dose of AT leaves the cardiac pacemaker under the sole sympathetic control, thus simplifying cardiac regulation and reducing the complexity of the cardiac control (Porta et al 2007c, Porta et al 2000). The β-adrenergic sympathetic blockade obtained via the administration of PR did not induce any variation on the complexity of cardiac regulation, thus suggesting that the contribution of sympathetic circuits to the complexity of cardiac control is negligible. The important reduction of the complexity of cardiac control after AT and the insignificant effect of PR were further confirmed by double-blockade of vagal and sympathetic branches obtained via the administration of PR after AT: indeed, CIs decreased after AT + PR and their reduction was similar to that after AT. The negligible contribution of sympathetic circuits to the complexity of cardiac control was further emphasized by the inability of CL to modify CIs and NCIs derived from the HP series.

The analysis of the complexity of the vascular control, as provided by CIs derived from the SAP series, sheds further light about cardiovascular regulation. The complexity of vascular control was not affected by high doses of AT (Porta et al 2000) as a likely result of the missing vagal influences over the vascular tree. In addition, the complexity of the vascular control was unmodified by PR, thus suggesting that sympathetic controls different from the β-adrenergic one are more effective in modulating SAP. This observation was confirmed by CL: indeed, the central blockade of all sympathetic influences affected the complexity of SAP. The increase of CIs after CL might be the effect of the desynchronizing effects over the vascular regulation induced by full sympathetic blockade. NCIs confirmed results obtained by CIs. It appears that normalization in this experimental protocol did not produce any additional advantage.

Despite the theoretical advantages of the KNNCE approach, CI_KNNCE (or NCI_KNNCE) showed the same ability as CI_CApEn and CI_SampEn (or NCI_CApEn and NCI_SampEn) to distinguish different experimental conditions. This result suggests that very few past samples were sufficient to predict future behaviors and their number remained constant over different experimental conditions. It is worth noting that KNNCE was more powerful in detecting the decrease in cardiac complexity after AT + PR than CCE. The limited discriminative power of CI_CCE and NCI_CCE might be the result of the larger dispersion of the indexes introduced by the corrective term.

The findings relevant to the KNNCE approach were confirmed when SE(y/y_L−1(i − 1)) was assessed as the negative natural logarithm of the probability that two samples of y/y_L−1(i − 1) are at distance closer than r = 20% of the standard deviation of y instead of computing SE(y/y_L−1(i − 1)) after UQ over ξ = 6 bins. This observation suggests that findings relevant to the KNNCE approach are independent of the strategy utilized to assess the SE of the conditional distributions.