Hypoglycemia-Induced Decrease of EEG Coherence in Patients with Type 1 Diabetes

Experimental protocol and database

The database is taken by a study approved by the local ethical committee, and data from this study have previously been published by Sejling et al. ⁶ and Fabris et al. ⁸ In particular, 19 patients with T1D participated in this study. Inclusion and exclusion criteria have previously been described by Sejling et al. ⁶ The diabetes patients underwent a hyperinsulinemic–hypoglycemic clamp with a target of 2.0–2.5 mmol/L. Hypoglycemia was induced by an intravenous infusion of fast-acting recombinant human insulin (Actrapid^®; Novo Nordisk, Bagsvaerd, Denmark). The insulin was administered intravenously at an infusion rate of 1 mU/kg/min, and target blood glucose (BG) levels were obtained by a variable infusion of glucose (200 mg/mL). BG level was determined every 5 min by a laboratory analyzer (model YSI 2300; Yellow Springs Instrument Co., Yellow Springs, OH).

EEG recordings and glucose concentration samples were collected in parallel during the experiment. Participants were sitting in an armchair, not allowed to sleep, and under constant observation by a physician who evaluated signs of neuroglycopenia, while 19 EEG channels were recorded by a digital EEG recorder (model Easy II; Cadwell Laboratories, Kennewick, WA) using standard cap electrodes placed on the scalp according to the 10/20 international system. In particular, Fp1-A1A2, Fp2-A1A2, F7-A1A2, F8-A1A2, F3-A1A2, F4-A1A2, Fz-A1A2, T3-A1A2, T4-A1A2, T5-A1A2, T6-A1A2, C3-A1A2, C4-A1A2, Cz-A1A2, P3-A1A2, P4-A1A2, Pz-A1A2, O1-A1A2, and O2-A1A2 EEG channels were acquired. (Odd numbers represent the left hemisphere, and even numbers represent the right hemisphere. The capital letters F, T, C, P, and O identify the frontal, temporal, central parietal, and occipital lobes respectively.) The EEG scan was analogically low-pass-filtered at 70 Hz and then digitally acquired and sampled at 200 Hz. The dynamic range of the EEG was ±4,620 μV with an amplitude resolution of 0.14 μV. Internal noise level in the analog data acquisition system was estimated to be 1.3 μV root mean square. Figure 1 shows the BG data (Fig. 1A), as well as the C4-A1A2 and O1-A1A2 EEG recordings (Fig. 1B and E) collected in parallel in a representative subject during the daytime.

Pre-analysis

Only EEG channels in the temporal, central, parietal, and occipital lobes (i.e., the L = 12 EEG time-series referred to, respectively, as the T3-A1A2, T4-A1A2, T5-A1A2, T6-A1A2, P3-A1A2, P4-A1A2, Pz-A1A2, C3-A1A2, C4-A1A2, Cz-A1A2, O1-A1A2, and O2-A1A2 channels) were considered because their analysis is simplified from the virtual absence of artifacts (e.g., eye-induced, muscle activation-induced). EEG data were not band-pass-filtered before computing iPDC in order to avoid possible distortions in the results. ⁹

A pre-analysis was developed considering, in each subject, two intervals referred to as eu- and hypoglycemia, respectively, targeted from the available BG data. In particular, for each subject, two 1-h intervals corresponding to eu- and hypoglycemia conditions, respectively, were identified from the BG time-series by visual inspection detecting when the glycemic thresholds at 70 and 180 mg/dL were crossed by a smoothing spline approximation of the samples (Fig. 1A). Representative randomly chosen 5-s EEG sweeps are also depicted to appreciate, albeit qualitatively, differences of the EEG signal in eu- (Fig. 1C and F) and hypoglycemia (Fig. 1D and G). Assessing the synchronization of pairs of EEG sweeps in the same glycemic condition by visual inspection is challenging, calling for the use of multivariate analysis tools.

iPDC definition

Coherence is a standard investigational tool in the multivariate analysis of biomedical signals, especially for EEG, and several definitions have been proposed. ^10,11 In particular, the iPDC has been recently presented by Takahashi et al., ¹² as an improvement of partial directed coherence ¹³ and general partial directed coherence. ¹⁴ Indeed, iPDC is considered more stable than other approaches in coherence assessment when signals have large differences in their amplitudes, as happens with EEG in hypoglycemia compared with EEG in euglycemia in our case study.

Roughly speaking, iPDC is a multivariate spectral measure to compute the directed influences between any given pair of signals (i,j) of a multivariate dataset. This information is condensed in a complex function \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$iPDC_{i \leftarrow j} ( f )$$ \end{document} of the frequency f, which measures the relative interaction of the signal j with regard to signal i as compared with interactions of all js with other signals in the multivariate dataset. Although we refer the reader to Takahashi et al. ¹² for the mathematical details, the procedure for computing iPDC is briefly described by the following two steps.

In the first step, an MAR model is identified by fitting it against EEG data streams acquired from L different channels. Formally, letting x ( n ) \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}%\begin{align*}%{\textbf{\textit{x}} ( \textbf{\textit{n}} ) } = [%{\textbf{\textit{x_1}} ( \textbf{\textit{n}} ) ,%\textbf{\textit{x_2}} ( \textbf{\textit{n}} ) , \ldots ,%\textbf{\textit{x_L}} ( \textbf{\textit{n}} ) } ] ^T \ \ n = 1 , .%, N \tag{1}%\end{align*} \end{document}

(where T stands for transpose) be the vector containing the samples at the n-th sampling instant n (n = 1, …N) of the L considered EEG channels x ₁ , x ₂ , …, x_L , the MAR model of order p is described as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{\textbf{\textit{x}} ( \textbf{\textit{n}} ) } = - { \Sigma}_{r = 1}^p \textbf{\textit{A}}_\textbf{\textit{rx}} ( \textbf{\textit{n- r}} ) + \bi w ( \bi n ) \tag{2} \end{align*} \end{document}

where A_r (r = 1, …, p) are the L × L unknown matrices of model coefficients, and the L-size column vector w ( n ) is the innovation process with covariance matrix Σ_W . ¹² The matrices A_r (r = 1, …, p) and Σ_W are estimated by least squares exploiting the so-called Yule–Walker equations ¹⁵ : \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\def\bi{\boldsymbol}\begin{document} \begin{align*} \begin{array}{ll} \left[ \begin{matrix} \textbf{\textit{R}}_{\textbf{\textit{{xx}}}} ( { \bf 0} ) \quad \textbf{\textit{R}}_{\textbf{\textit{{xx}}}} ( { \bf - 1} ) & \ldots & \textbf{\textit{R}}_{\textbf{\textit{{xx}}}} ( \textbf{\textit{- p + {\bf 1}}} ) \\ \textbf{\textit{R}}_{\textbf{\textit{{xx}}}} ( { \bf 1} ) \quad \textbf{\textit{R}}_{\textbf{\textit{{xx}}}} ( { \bf 0} ) & & \textbf{\textit{R}}_{\textbf{\textit{{xx}}}}( \textbf{\textit{- p + {\bf 2}}} ) \\ \vdots & \ddots & \vdots \\ \textbf{\textit{R}}_\textbf{\textit{{xx}}} ( \textbf{\textit{p}} -{\bf 1} ) \quad \textbf{\textit{R}}_\textbf{\textit{{xx}}} ( \textbf{\textit{p}} - {\bf 2} ) & \cdots & \textbf{\textit{R}}_{\textbf{\textit{{xx}}}} ( { \bf 0} ) \end{matrix} \right]&\left[ \begin{matrix} \textbf{\textit{A}}^{\textbf{\textit{T}}}_{\bf 1} \\ \textbf{\textit{A}}^{\textbf{\textit{T}}}_{\bf 2} \\ \vdots \\ \textbf{\textit{{A}}}^{\textbf{\textit{T}}}_{\textbf{\textit{p}}} \end{matrix} \right] \\ \quad= - \left[ \begin{matrix} \textbf{\textit{R}}_{\textbf{\textit{{xx}}}} ( { \bf 1} ) \\ \textbf{\textit{R}}_{\textbf{\textit{xx}}} ( { \bf 2} ) \\ \vdots \\ \textbf{\textit{R}}_{\textbf{\textit{xx}}} ( \textbf{\textit{p}} ) \end{matrix} \right]& \tag{3} \end{array} \end{align*} \end{document}

where the R _xx ( k ) matrices are the covariance matrices computed from the available EEG data as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \textbf{\textit{R}}_{\textbf{\textit{{ xx } } }} ( \textbf { \textit { k } } ) = \begin{cases} \frac { 1 } { N } \Sigma_ { n = 0 }^{N - 1 - k } \textbf{\textit{ x } } ( \textbf{\textit { n } } ) \textbf{\textit{ x}}^{\textbf{\textit{T } }} ( \textbf{\textit { n + k} } ) \quad \quad \quad \quad \quad \quad k = 0 , 1 , . , N - 1 \\ \frac{1} { N } \Sigma_ { n = - k } ^ { N - 1 } \textbf{\textit{x} } ( \textbf{\textit{n}} ) \textbf{\textit{x}}^{\textbf{\textit{T} }} ( \textbf{\textit{n + k } } ) \ \ k = - ( N - 1 ) , - ( N - 2) , . , - 1 \end{cases}\tag{4} \end{align*} \end{document}

In Eq. 2 model order p is found by trials until the Akaike information criterion (AIC) index \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}AIC ( p ) = Nln \mid { \boldsymbol{\Sigma}}_\textbf{\textit{W}} \mid + 2L^2 p \tag{5} \end{align*} \end{document}

(where | | denotes the matrix determinant) is minimized. ¹⁶

In the second step, having defined the complex matrix B ( f ) as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \textbf{\textit{B}} ( \textbf{\textit{f}} ) = \textbf{\textit{I}}_{\textbf{\textit{L}}} - \Sigma_{r = 1}^p \textbf{\textit{A}}_{\textbf{\textit{r}}} e^{ - j2 \pi f} \tag{6} \end{align*} \end{document}

where I_L is the identity matrix and j is the imaginary unit in this equation, the iPDC complex function from the time-series j to the time-series i is obtained by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} iPDC_{i \leftarrow j} (f)= \sigma_{W_{ii}}^{-1/2} \frac{b_ij(f)}{\sqrt{\textbf{\textit{b}}_{\textbf{\textit{j}}}^{\textbf{\textit{H}}}( f){\boldsymbol{\Sigma}}_{\textbf{\textit{W}}}^{\textbf{\textit{- 1}}} \textbf{\textit{b}}_{\textbf{\textit{j}}} (f)}} \tag{7} \end{align*} \end{document}

where b_j (f) and b_ij(f) are, respectively, the j-th column and the (j,i)-th element of matrix B ( f ), \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\sigma_{W_{ii}}$$ \end{document} is the (i,i)-th element of the innovation covariance matrix Σ_W , and the apex H in \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\textbf{\textit{b}}_{\textbf{\textit{j}}}^{\textbf{\textit H}}$$ \end{document} stands for the Hermitian transpose (i.e., obtained from b_j by taking the transpose and then the complex conjugate of its components).

The complex function \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$iPDC_{i \leftarrow j} ( f )$$ \end{document} of Eq. 7 is usually analyzed in terms of its absolute value. ¹²

iPDC computation

iPDC analysis was applied to the first 15 min of the eu- and hypoglycemia intervals. As previously described in the Pre-analysis section, L = 12 EEG channels were considered, and each 15-min interval of eu- and hypoglycemia was divided in 5-s sweeps, resulting in N = 1,000 samples per sweep, a quantity sufficient to obtain stable results in terms of coherence estimation according to Florin et al. ⁹ At the same time, 5-s intervals are narrow enough to invoke pseudostationarity assumptions on EEG (e.g., as done by Abásolo et al. ¹⁷ ). For each sweep, an MAR model as in Eq. 2 was identified. Regarding the model order p, we first considered, for each sweep, the value minimizing AIC(p) of Eq. 5 among candidate orders ranging from 1 to 20; then, among the estimated values of p, we chose the highest and used it as the model order for all sweeps under analysis for a given subject (the resulting p was always 9 or 10). The rationale behind this decision is that, to detect coupling directions in multivariate oscillatory systems, the MAR model should have a high order, ¹⁸ and at the same time, overestimating the model order is better than underestimating it. ⁹ Then, for any pair (i,j), the function \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$iPDC_{i \leftarrow j} ( f )$$ \end{document} of Eq. 7 was computed in theta ([4–8] Hz) and alpha ([8–13] Hz) bands for each sweep. Eventually, the absolute mean profiles \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\overline{PDC}}_{i \leftarrow j} ( f )$$ \end{document} , referred to the whole interval of 15 min, were computed by averaging the absolute values of the \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$iPDC_{i \leftarrow j} ( f )$$ \end{document} profiles estimated in the 5-s sweeps (throughout this article, for sake of brevity, the symbol \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\overline{{PDC}}}$$ \end{document} will be used without explicitly showing the dependence on frequency and channel pair).

iPDC values analysis

Because the values of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \overline {{PDC}}}$$ \end{document} can be pretty low, it is necessary, before evaluating changes passing from eu- to hypoglycemia, to determine an analytical threshold exploiting surrogate data. In particular, 15-min surrogate data in eu- and hypoglycemia were obtained randomly changing the time order of a set of real data (i.e., randomly permuting in temporal order the samples of the original series and so maintaining mean, variance, and histogram distribution of the original data). ¹⁹ Then, as for real data, iPDC was estimated, and two thresholds for eu- and hypoglycemia, respectively, were calculated as 5% of the absolute value of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \overline {{PDC}}}$$ \end{document} . ²⁰ Thus, before applying the statistical tests to evaluate if there were remarkable changes passing from eu- to hypoglycemia, it was checked if \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \overline {{PDC}}}$$ \end{document} results were above the threshold estimated at any frequency. Eventually, to assess for each subject whether or not the changes in the absolute value of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\overline{PDC}}$$ \end{document} were statistically significant when passing from eu- to hypoglycemia, paired Student's t tests were computed under the hypothesis of normal distribution of samples (Lilliefors test); otherwise, Wilcoxon rank-sign tests were considered (a value of P < 0.05 was considered significant). After that, for each subject, the average value of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\overline{PDC}}$$ \end{document} was computed to obtain a scalar indicator for the eu- and hypoglycemia interval. This scalar indicator (which subsequently in this article will be denoted by the symbol \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${\overline{\overline {{{PDC}}}}}$$ \end{document} ) is obviously much simpler to handle than the function \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$${ \overline{PDC}}$$ \end{document} for the analysis.