A probabilistic framework for a physiological representation of dynamically evolving sleep state

Abstract

This work presents a probabilistic method for mapping human sleep electroencephalogram (EEG) signals onto a state space based on a biologically plausible mathematical model of the cortex. From a noninvasive EEG signal, this method produces physiologically meaningful pathways of the cortical state over a night of sleep. We propose ways in which these pathways offer insights into sleep-related conditions, functions, and complex pathologies. To address explicitly the noisiness of the EEG signal and the stochastic nature of the mathematical model, we use a probabilistic Bayesian framework to map each EEG epoch to a distribution of likelihoods over all model sleep states. We show that the mapping produced from human data robustly separates rapid eye movement sleep (REM) from slow wave sleep (SWS). A Hidden Markov Model (HMM) is incorporated to improve the path results using the prior knowledge that cortical physiology has temporal continuity.

Appendix A: Mathematical model of sleep

The dimensionless version of the model detailed in Kramer et al. (2005) for seizures and in Lopour et al. (2011) for sleep. In Tables 1 and 2, the definitions of the variables and parameters of the model are described.

Table 1 Model variables

Table 2 Model Parameters

$$ \frac{\partial\tilde{h}_{e}}{\partial\tilde{t}}=1-\tilde{h}_e+ \frac{\Delta h_{e}^{rest}}{h^{rest}}+ L\Gamma_{e}(h_{e}^{0}-\tilde{h_{e}})\tilde{I}_{ee} +\Gamma_{i}(h_{i}^{0}-\tilde{h_{e}})\tilde{I}_{ie} $$

(7)

$$ \frac{\partial\tilde{h}_{i}}{\partial\tilde{t}}=1-\tilde{h}_i+ L\Gamma_{e}(h_{e}^{0}-\tilde{h_{i}})\tilde{I}_{ei} +\Gamma_{i}(h_{i}^{0}-\tilde{h_{i}})\tilde{I}_{ii} $$

(8)

$$ \left(\frac{1}{T_{e}}\frac{\partial}{\partial\tilde{t}}+1\right)^{2} \tilde{I}_{ee}= N_{e}^{\beta}\tilde{S}_{e}\left[\tilde{h}_{e}\right] +\tilde{\phi}_e+P_{ee}+\tilde{\Gamma}_1 $$

(9)

$$ \left(\frac{1}{T_{e}}\frac{\partial}{\partial\tilde{t}}+1\right)^{2} \tilde{I}_{ei}= N_{e}^{\beta}\tilde{S}_{e}\left[\tilde{h}_{e}\right] +\tilde{\phi}_i+P_{ei}+\tilde{\Gamma}_2 $$

(10)

$$ \left(\frac{1}{T_{i}}\frac{\partial}{\partial\tilde{t}}+1\right)^{2} \tilde{I}_{ie}= N_{i}^{\beta}\tilde{S}_{i}\left[\tilde{h}_{i}\right] +P_{ie}+\tilde{\Gamma}_3 $$

(11)

$$ \left(\frac{1}{T_{i}}\frac{\partial}{\partial\tilde{t}}+1\right)^{2} \tilde{I}_{ii}= N_{i}^{\beta}\tilde{S}_{i}\left[\tilde{h}_{i}\right] +P_{ii}+\tilde{\Gamma}_{4} $$

(12)

$$ \left(\frac{1}{\lambda_{e}}\frac{\partial}{\partial\tilde{t}}+1\right)^{2} \tilde{\phi}_{e}= \frac{1}{\lambda_{e}^{2}} \frac{\partial^{2}\tilde{\phi}_{e}}{\partial\tilde{x}^{2}}+ \left(\frac{1}{\lambda_{e}}\frac{\partial}{\partial \tilde{t}}+1\right)N_{e}^{\alpha}\tilde{S}_{e}\left[\tilde{h}_{e}\right] $$

(13)

$$ \left(\frac{1}{\lambda_{i}}\frac{\partial}{\partial\tilde{t}}+1\right)^{2} \tilde{\phi}_{i}= \frac{1}{\lambda_{i}^{2}}\frac{\partial^{2}\tilde{\phi}_{i}}{\partial\tilde{x}^{2}}+ \left(\frac{1}{\lambda_{i}}\frac{\partial}{\partial \tilde{t}}+1\right)N_{i}^{\alpha}\tilde{S}_{e}\left[\tilde{h}_{e}\right] $$

(14)

The sigmoid transfer functions that transform the mean soma voltage to the dimensionless firing rates ($\tilde {S}_{e}, \ \tilde {S}_{i}$) are shown in Eqs. (15) and (16).

$$ \tilde{S}_{e}\left[\tilde{h}_{e}\right] = \frac{1}{1+\exp\left[-\tilde{g}_{e}(\tilde{h}_e-\tilde{\theta}_e)\right]} $$

(15)

$$ \tilde{S}_{i}\left[\tilde{h}_{i}\right] = \frac{1}{1+\exp\left[-\tilde{g}_{i}(\tilde{h}_i-\tilde{\theta}_i)\right]} $$

(16)

The stochastic input from subcortical sources is:

$$\begin{array}{@{}rcl@{}} \tilde{\Gamma}_1&=&\alpha_{ee}\sqrt{P_{ee}}\xi_1[\tilde{x},\tilde{t}] \nonumber\\ \tilde{\Gamma}_2&=&\alpha_{ei}\sqrt{P_{ei}}\xi_2[\tilde{x},\tilde{t}] \nonumber \\ \tilde{\Gamma}_3&=&\alpha_{ie}\sqrt{P_{ie}}\xi_3[\tilde{x},\tilde{t}] \nonumber \\ \tilde{\Gamma}_4&=&\alpha_{ii}\sqrt{P_{ii}}\xi_4[\tilde{x},\tilde{t}] \nonumber\\ \alpha_{ee, ei, ie, ii}&=&\hat{\alpha}(-4\mu_{e} + 5) \end{array} $$

(17)

where ξ _k are white Gaussian noise with zero mean, δ-function correlations and approximated numerically by:

$$ \xi_{k}[\tilde{x},\tilde{t}] = \frac{R(m,n)}{\sqrt{\Delta\tilde{x}\Delta{\tilde{t}}}} $$

(18)

and $[\tilde {x},\tilde {t}] = [m\Delta {\tilde {x}},\ \ n\Delta {\tilde {t}}]$–integer coordinates on a grid spaced by $[\Delta {\tilde {x}},\ \ \Delta {\tilde {t}}]$. Finally, Eq. (17) for α _{ee, ei, ie, ii} is from Lopour et al. (2011). This PDE was simulated in MATLAB with an Euler predictor-corrector method designed in Lopour et al. (2011).

Appendix B: Features

B.1 Initial feature set

• Power in the delta, theta, alpha, beta, and gamma frequency ranges, and total power

• Statistical properties: variance, skewness, and kurtosis

• Composite permutation entropy index (CPEI) (Olofsen et al. 2008)

• Spindle index (MacKay et al. 2010)

• Power fractions: high power fraction (beta and gamma fraction of total power), and low power fraction (delta and theta fraction of total power)

• Properties of the log power spectrum curve: slope and offset of a fitted line, maximum alpha power above the linear fit, and maximum delta power above the linear fit (Leslie et al. 2009; Lopour et al. 2011)

• Delta wave steepness properties: the average steepness of delta waves, fraction of epoch covered by delta waves, and combination delta wave steepness. The code for calculating this was developed for this paper by the authors and is included in Online Resource 2, a zip file containing the m-files deltaWaveSteepnessFxn.m and epochFeatsFromWaveSlope.m. However, the importance and meaning of this feature is highlighted in Riedner et al. (2007).

• Median frequency (Tsimpiris and Kugiumtzis 2010)

• Autocorrelation functions with delay τ = 1 time-step: Spearman autocorrelation, Pearson autocorrelation, and partial autocorrelation (Tsimpiris and Kugiumtzis 2010)

• Bicorrelation with delay τ = 1 time-step (Tsimpiris and Kugiumtzis 2010)

• Mutual information measures including: equidistant mutual information, equiprobable mutual information, and the first minimums of both types of mutual information. Algorithms were from (Tsimpiris and Kugiumtzis 2010), and for these features, with delay τ = 1 time-step

• Hurst exponent (Tsimpiris and Kugiumtzis 2010)

• Hjorth parameters (Tsimpiris and Kugiumtzis 2010)

Note that several of the more complex feature calculation algorithms were implemented using MATLAB code created by other researchers. In particular, many standard time series analysis features were taken from the MATS toolbox (Tsimpiris and Kugiumtzis 2010), and others from supplemental material to papers including Olofsen et al. (2008), MacKay et al. (2010), Leslie et al. (2009), Lopour et al. (2011). The code sources are indicated in the above listing.

B.2 Feature selection

The requirements stated in Section 4 provide a baseline for creating heuristic scoring functions to prune features and then select a limited number of feature subsets that are likely to perform well. First, the staged patient data is used to split epochs into SWS or REM sleep. The feature selection algorithm is the only time we make use of labeled sleep stage data in this algorithm other than for visualization purposes. Because both the model and patient data contain labels of areas or epochs that are SWS and REM sleep, this is applicable to both types of data and between the two types of data.

Then, with these labeled epochs, the Spearman rank correlation between SWS and REM of each feature was calculated in both the model and the sleep data. That is to say that SWS and REM were converted into a binary variable, respectively 0 and 1. Then the Spearman rank correlation was calculated between this binary variable and the feature values at a particular feature, giving both a correlation and a p-value. One can think of this as a measure of how likely an REM epoch feature value is to be greater than a SWS epoch feature value. If this correlation is close to 1, that means that most REM feature values are larger than SWS feature values. Conversely, a correlation close to -1 means that most SWS feature values are larger than REM feature values. A correlation close to 0 indicates that SWS and REM features values are well mixed in ranking–a typical feature value of a SWS epoch is equally likely to be greater than or less than a typical feature value of an REM epoch. If all SWS features are greater than all REM features, the correlation would be precisely negative one. The Spearman correlation is a non-linear measure of correlation, and only compares the ranking of values.

To remove features that did not vary significantly between SWS and REM in the patient data, all features with less absolute value Spearman correlation than a threshold of 0.35 in all patients were removed. The features removed due to lack of importance were the low power fraction, the skewness, gamma power, median frequency, bicorrelation, and the mutual information first minimums.

Next, using the same Spearman correlations, all remaining features that had opposite signs of correlation coefficient between the model and sleep were removed from the feature set if absolute value of the model Spearman correlation was greater than .25. This was implemented because a different sign of correlation coefficient is an indication that the model features vary with an opposite trend when compared to patient data. Clearly, opposite trends between SWS and REM demonstrate that the feature is not modeled well with the model, but only if the model shows a clear trend is this significant. Only few features were of opposite sign, and only the power in the β-band and the maximum alpha band power spectral density above a linear approximation exceeded that threshold.

Finally, to choose subsets of features that reduce redundancy, the Pearson correlation (a linear correlation) of features across the model for each state was calculated. Features that highly correlated to one another in the model space were identified by thresholding correlations and identifying features that had over ρ = . 9 correlation. From this, a connectivity graph of which features were connected by high correlations could be drawn, for example see Fig. 12 which shows one of the complete subgraphs, in which all features are redundant with each other and thus all features are fully connected to each other. To remove excess redundancy, the feature selection algorithm allowed only one feature from each complete subgraph to be contained in the final feature set. This also reduced the number of possible subsets.

Naive Bayes models also assume conditional independence of features (conditioned on the state), a ‘naive’ assumption that gives this model its name. We do not expect the Naive Bayes model assumptions to be met exactly, but to avoid high correlations, we do not allow feature sets to contain any two features that given the state have an average Pearson correlation coefficient of 0.4.

A test set of all feature subsets from size one to six were generated after pruning features and limiting redundancy according to these rules.

All subsets generated using the rules above were tested against a cost function designed to separate REM and SWS regions. As a significant goal is to distinguish SWS and REM when the EEG epochs are mapped to the sleep parameter space, the average likelihood of each stage of each subject for each feature subset was calculated. Figure 13 shows an example of this type of region.

Next, the states in parameter space that contained 95 % of the likelihood for average REM likelihood and SWS sleep were each identified. Areas to the left (lower L) and below (lower $\Delta h_{e}^{rest}$) the cusp of the manifold region were identified as SWS regions, while REM regions were identified as regions above the manifold region (higher $\Delta h_{e}^{rest}$). The manifold itself and the top left corner beyond the cusp of the manifold were allowed to be either REM or SWS - this is because the manifold is thought to be a transition region, and while low L and high $\Delta h_{e}^{rest}$ is theorized to contain REM, the results in Lopour et al. (2011) have indicated that perhaps some SWS may lie in these areas.

From these 95 % likelihood regions (composed of a set of states), the percentage overlap of REM into SWS-identified areas and SWS into REM-identified areas of the parameter space were calculated. A cost function composed of the sum of these percentages was minimized at each feature set size. Under visual inspection, the solution for the feature set with four features was chosen.