Baroreflex analysis in diabetes mellitus: linear and nonlinear approaches

Abstract

The aim of our study was to employ novel nonlinear synchronization approaches as a tool to detect baroreflex impairment in young patients with subclinical autonomic dysfunction in Type 1 diabetes mellitus (DM) and compare them to standard linear baroreflex sensitivity (BRS) methods. We recorded beat-to-beat pulse interval (PI) and systolic blood pressure (SBP) in 14 DM patients and 14 matched healthy controls. We computed the information domain synchronization index (IDSI), cross-multiscale entropy, joint symbolic dynamics, information-based similarity index (IBSI) in addition to time domain and spectral measures of BRS. This multi parametric analysis showed that baroreflex gain is well-preserved, but the time delay within the baroreflex loop is significantly increased in patients with DM. Further, the level of similarity between blood pressure and heart rate fluctuations was significantly reduced in DM. In conclusion, baroreflex function in young DM patients is changed. The quantification of nonlinear similarity and baroreflex delay in addition to baroreflex gain may provide an improved diagnostic tool for detection of subclinical autonomic dysfunction in DM.

Appendix

1.1 Information domain synchronization index (IDSI)

The first step in the algorithm is the normalization (subtracting the mean and dividing by the standard deviation) of both signals and their symbolization into symbols. We used six quantization levels (see “Methods” section).

Given a pair of normalized symbolized signals (x, y) with x = {x(i), i = 1,…, N} and y = {y(i), i = 1,…, N}, the cross-conditional entropy of x given a pattern of $y\left( {{\text{CE}}_{x/y} } \right)$ is defined as

$$ {\text{CE}}_{x/y} \left( L \right) = - \sum\limits_{L - 1} {p\left( {y_{L - 1} } \right)\sum\limits_{i/L - 1} {p\left( {x(i)/y_{L - 1} } \right){ \log }\,p\left( {x(i)/y_{L - 1} } \right)} } $$

where log performs the natural logarithm, p(y _L−1) denotes the probability of the pattern y _L−1 and $p\left( {x(i)/y_{L - 1} } \right)$ symbolizes the conditional probability of the sample x(i) given the pattern y _L−1. This represents the amount of information carried by the present sample of the signal x (i.e., x(i)) when a pattern of L − 1 samples of the signal y (i.e., $y_{L - 1} \left( i \right) = \left( {y(i), \ldots ,y(i - L + 2)} \right)$) is known. Similarly, ${\text{CE}}_{y/x} \left( L \right)$ is calculated for various lengths of pattern (in our case L varied from 1 to 12).

In the next step, the normalized cross-conditional entropy ${\text{NCE}}_{x/y} \left( L \right)$ and ${\text{NCE}}_{y/x} \left( L \right)$ are obtained by normalizing ${\text{CE}}_{x/y} \left( L \right)$ and ${\text{CE}}_{y/x} \left( L \right)$ by E _x (1) and E _y (1), respectively. E _x (1) and E _y (1)denote the Shannon entropies of x and y, respectively

$$ E_{x} (1) = - \sum\limits_{L} {p(x(i)){ \log }\,p(x(i))} $$

$$ E_{y} (1) = - \sum\limits_{L} {p(y(i)){ \log }\,p(y(i))} $$

representing the amount of information carried by the sample x(i) (or y(i)) when no previous sample is given.

Since both ${\text{NCE}}_{x/y} \left( L \right)$ and ${\text{NCE}}_{y/x} \left( L \right)$ decrease towards zero as an effect of the shortness of the data, they are substituted by normalized corrected cross-conditional entropies ${\text{NCCE}}_{x/y} \left( L \right)$ and ${\text{NCCE}}_{y/x} \left( L \right)$ calculated as

$$ {\text{NCCE}}_{x/y} \left( L \right) = {\text{NCE}}_{x/y} \left( L \right) + {\text{perc}}_{y} \left( L \right) $$

$$ {\text{NCCE}}_{y/x} \left( L \right) = {\text{NCE}}_{y/x} \left( L \right) + {\text{perc}}_{x} \left( L \right) $$

where perc _y (L) and perc _x (L) are percentages of length L patterns found only one time in the data set y and x, respectively.

Finally, as a measure of two signals uncoupling the uncoupling function (UF) is calculated

$$ {\text{UF}}_{x,y} \left( L \right) = { \min }\left( {{\text{NCCE}}_{y/x} \left( L \right),{\text{NCCE}}_{x/y} \left( L \right)} \right) $$

where min takes the minimum value between ${\text{NCCE}}_{y/x}$ and ${\text{NCCE}}_{x/y}$ for each L. The larger the index is the more independent is the signals. An index of synchronization (IDSI) quantifies the maximum amount of information exchange between two signals and is defined as

$$ {\text{IDSI}} = 1 - { \min }\left( {{\text{UF}}_{x/y} \left( L \right)} \right) $$

1.2 Cross multiscale entropy (Cross-MSE)

Cross multiscale entropy analysis is the analysis of asynchrony of two signals for various time scales; it is based on the calculation of Cross Sample Entropy (Cross-SampEn) values as a function of time scale.

For computation of Cross Sample Entropy (Cross-SampEn), three input parameters m, r, and N have to be fixed. From two normalized (subtracting the mean and dividing by the standard deviation) time series x and y (x = {x(i), i = 1,…, N} and y = {y(i), i = 1,…, N}), vectors x _m (i) and y _m (i) are formed, where x _m (i) = {x(i + k): 0 ≤ k ≤ m − 1} and y _m (i) = {y(i + k): 0 ≤ k ≤ m − 1} are the vectors of m consecutive data points of corresponding signal.

The distance between two such vectors is defined to be d[x _m (i), y _m (j)] = max {| x(i + k) − y(j + k) |: 0 ≤ k ≤ m − 1}, i.e., the maximum difference of their corresponding scalar components.

The whole set of vector sequences x _m (1), x _m (2),…, x _m (N − m) is stepwise used to calculate values $B_{i}^{m} (r)(y\parallel x)$ as (N − m)⁻¹ times the number of vectors y _m (j) within distance r of x _m (i), where j ranges from 1 to N − m. Then, B ^m(r) is defined as:

$$ B^{m} (r)(y\parallel x) = (N - m)^{ - 1} \sum\limits_{i = 1}^{N - m} {B_{i}^{m} } (r)(y\parallel x) $$

Similarly, $A_{i}^{m} (r)(y\parallel x)$ is defined as (N − m)⁻¹ times the number of vectors y _m+1(j) within r of x _m+1(i), where j ranges from 1 to N − m, and set

$$ A^{m} (r)(y\left\| {x)} \right. = (N - m)^{ - 1} \sum\limits_{i = 1}^{N - m} {A_{i}^{m} } (r)(y\left\| {x)} \right. $$

The Cross-SampEn(m, r, N) is then defined as

$$ {\text{CrossSampEn}}(m,r,N) = - { \ln }\left[ {A^{m} (r)/B^{m} (r)} \right] $$

Cross-multiscale entropy analysis was computed according to the concept proposed by Costa [9]. Given a one-dimensional discrete time series, {x ₁,…, x _i ,…, x _N }, we constructed consecutive coarse-grained time series {x ^(τ)} determined by the scale factor τ, according to the equation:

$$ x_{j}^{(\tau )} = 1/\tau \sum\limits_{i = (j - 1)\tau + 1}^{j\tau } {x_{i} } $$

where τ represents the scale factor and 1 ≤ j ≤ N/τ. In other words, coarse-grained time series for scale τ were obtained by taking arithmetic mean of τ neighboring original values without overlapping. Similarly, {y ^(τ)} for signal y were calculated. For scale 1, the coarse-grained time series is simply the original time series. The CrossSampEn values are calculated for each of the coarse-grained pairs (for given τ) of signals.

1.3 Joint symbolic dynamics

In X (Eq. 1), x ^PIand x ^SBP are n beat-to-beat values of PI and SBP, respectively:

$$ {\bf X} = \left\{ {\left[ {x_{n}^{\text{PI}} ,x_{n}^{\text{SBP}} } \right]^{T} } \right\}_{n = 0,1, \ldots } x \in R $$

(1)

X is transformed in S defined as:

$$ {\bf S} = \left\{ {\left[ {s_{n}^{\text{PI}} ,s_{n}^{\text{SBP}} } \right]^{T} } \right\}_{n = 0,1,..} s \in 0,1 $$

with the following definitions:

$$ s_{n}^{\text{PI}} = \left\{ \begin{gathered} 0:(x_{n}^{\text{PI}} - x_{n + 1}^{\text{PI}} ) \le l^{\text{PI}} \\ 1:(x_{n}^{\text{PI}} - x_{n + 1}^{\text{PI}} ) > l^{\text{PI}} \\ \end{gathered} \right.\quad \quad s_{n}^{\text{SBP}} = \left\{ \begin{gathered} 0:(x_{n}^{\text{SBP}} - x_{n + 1}^{\text{SBP}} ) \le l^{\text{SBP}} \\ 1:(x_{n}^{\text{SBP}} - x_{n + 1}^{\text{SBP}} ) > l^{\text{SBP}} \\ \end{gathered} \right. $$

where threshold value l is set zero. Thus, increases between two successive PI and SBP, respectively are coded as ‘1’ and consequently decreases and equilibrium are coded as ‘0’. Subsequently, S is subdivided into short sequences with a certain length. Each single word is obtained by a shift of one within the symbol string S. The length of words is limited due to the requirement of a statistically sufficient representation of each single word type. For 60-min recordings there are no more than 64 different word types feasible. Taking into account the four different symbol combinations within S (the alphabets of PI as well as SBP consist each of two elements), words with a maximum length of three are possible (2³*2³ = 64). Therefore, this approach is able to map the dynamics of PI and SBP within four consecutive heart beats (i.e., three PI). From the perspective of phase space, this corresponds to a three-dimensional embedding, which is a pragmatic rather then a true reconstruction of the system.

The word distribution density matrix W contains the frequency of each of the 8 × 8 possible combinations of PI and SBP patterns.

$$ {\bf W} = \left[ {\begin{array}{*{20}c} {PI_{000} ,SBP_{000} } & \cdots & {PI_{000} ,SBP_{111} } \\ \vdots & \ddots & \vdots \\ {PI_{111} ,SBP_{000} } & \cdots & {PI_{111} ,SBP_{111} } \\ \end{array} } \right] $$

Statistical analysis [1] showed that some word types within W are significantly over-represented and other ones are under-represented. All word types representing symmetric PI and SBP behavior (i.e., W _1,1, W _2,2, W _3,3, W _4,4, W _5,5, W _6,6, W _7,7, W _8,8) occur frequently, whereas some word types representing diametric behavior (i.e. W _7,2, W _5,4, W _4,5, W _2,7) occur significantly less often. Interpreting the former ones as baroreflex patterns and latter ones as missing or a lack of baroreflex responses, the diagonals of the matrix were summed up: ${\text{JSDsym}} = \frac{1}{N}\sum\nolimits_{j = k = 1}^{8} {W_{j,k} }$—representing symmetric word types and ${\text{JSDdiam}} = \frac{1}{N}\sum\nolimits_{j = k = 1}^{8} {W_{j,9 - k} }$—representing diametric word types within W.

1.4 Information-based similarity index

The first step of computing the IBSI involves symbolization. Changes in successive values of x and y, respectively, are encoded as ‘1’ if an increase and as ‘0’ if a decrease in respective time series occurred. We considered words of k-bit length, shifting one data point at a time. We counted the occurrences of different word types, and subsequently sorted them in descending order by frequency of occurrence. The resulting rank-frequency distribution represents the statistical hierarchy of symbolic words. Therefore, the first rank word corresponds to the type of fluctuation which is the most frequent pattern in the time series.

In the next step, the rank order difference between x and y time series was visualized by plotting the rank number of each k-bit word in the first time series against that of the second time series. The dissimilarity between x and y time series was quantified by measuring the scatter of these points. If two time series are similar in their rank order of the words, the scattered points will be located near the diagonal line. Therefore, the average deviation of these scattered points away from the diagonal line can provide information about the similarity between two series and is expressed as IBSI value within the interval [0; 1]. IBSI is calculated as a weighted deviation between two symbolic sequences x and y:

$$ {\text{IBSI}} = {\frac{1}{{2^{m} - 1}}}\sum\limits_{k = 1}^{{2^{m} }} {\left| {R_{x} (w_{k} ) - R_{y} (w_{k} )} \right|} F(w_{k} ) $$

where

$$ F(w_{k} ) = \frac{1}{Z}\left[ { - p_{x} \left( {w_{k} } \right){ \log }\,p_{x} \left( {w_{k} } \right) - p_{y} \left( {w_{k} } \right){ \log }\,p_{y} \left( {w_{k} } \right)} \right] $$

Here p _x (w _k ) and R _x (w _k ) represent probability and rank of a specific word w _k , in time series x. Similarly, p _y (w _k ) and R _y (w _k ) stand for probability and rank of the same k-bit word in time series y. The normalization factor Z is given by

$$ Z = \sum\nolimits_{k} {[ - p}_{x} \left( {w_{k} } \right){ \log }\,p_{x} \left( {w_{k} } \right) - p_{y} \left( {w_{k} } \right){ \log }\,p_{y} \left( {w_{k} } \right)] $$