Modeling the Complex Dynamics and Changing Correlations of Epileptic Events

To scale our model to a large number of channels, we consider a Gaussian graphical model (GGM) for $\boldsymbol{\rm\epsilon}_{t}$ capturing a sparse dependency structure amongst the channels. Let $G=(V,E)$ be an undirected graph with $V$ the set of channel nodes $i$ and $E$ the set of edges with $(i,j)\in E$ if $i$ and $j$ are connected by an edge in the graph. Then, $\left[\Delta_{l}^{-1}\right]_{ij}=0$ for all $(i,j)\not\in E$, implying $\epsilon_{t}^{(i)}$ is conditionally independent of $\epsilon_{t}^{(j)}$ given $\epsilon_{t}^{(k)}$ for all channels $k\neq i,j$. In our dynamic model of Eq. (1), statements of conditional independence of $\boldsymbol{\rm\epsilon}_{t}$ translate directly to statements of the observations $\boldsymbol{\rm y}_{t}$.

In our application, we choose $G$ based on the spatial adjacencies of channels in the electrode grid, as depicted in Figure 1 (bottom left). In addition to encoding the spatial proximities of iEEG electrodes, the graphical model also yields a sparse precision matrix $\Delta_{l}^{-1}$, allowing for more efficient scaling to the large number of channels commonly present in iEEG. These computational efficiencies are made clear in Section 3.

Recall that our formulation involves $N+1$ independently evolving Markov chains: $N$ chains for the channel states $z_{t}^{(i)}$ plus one for the event state sequence $Z_{t}$. As indicated by the observation model of Eq. (3), the $N+1$ Markov chains jointly generate our observation vector $\boldsymbol{\rm y}_{t}$ leading to an interpretation of our formulation as a factorial HMM [14]. However, here we have a sparse dependency structure in how the Markov chains influence a given observation $\boldsymbol{\rm y}_{t}$, as induced by the conditional independencies in $\boldsymbol{\rm\epsilon}_{t}$ encoded in the graph $G$. That is, $y_{t}^{(i)}$ only depends on $Z_{t}$ the set of $z_{t}^{(j)}$ for which $j$ is a neighbor of $i$ in $G$.

As in the AR-HMM, we place a multivariate normal prior on the AR coefficients,

with mean $\boldsymbol{\rm m}_{0}$ and covariance $\Sigma_{0}$. Throughout this work, we let $\boldsymbol{\rm m}_{0}=\boldsymbol{\rm 0}$.

For the channel covariances $\Delta_{l}$ with sparse precisions $\Delta_{l}^{-1}$ determined by the graph $G$, we specify a hyper-inverse Wishart (HIW) prior,

where $b_{0}$ denotes the degrees of freedom and $D_{0}$ the scale. The HIW prior [15] enforces hyper-Markov conditions specified by $G$.

A natural question is how many AR states are the channels switching between? Likewise, which are shared between the channels and which are unique? We expect to see similar dynamics present in the channels (sharing of AR processes), but also some differences. For example, maybe only some of the channels ever get excited into a certain state. To capture this structure, we take a Bayesian nonparametric approach building on the beta process (BP) AR-HMM of Fox et al. [16]. Through the beta process prior [17], the BP-AR-HMM defines a shared library of infinitely many AR coefficients $\{\boldsymbol{\rm a}_{k}\}$, but encourages each channel to only use a sparse subset of them.

The BP-AR-HMM specifically defines a featural model. Let $\boldsymbol{\rm f}^{(i)}$ be a binary feature vector associated with channel $i$ with $f_{k}^{(i)}=1$ indicating that channel $i$ uses the dynamic described by $\boldsymbol{\rm a}_{k}$. Formally, the feature assignments $f_{k}^{(i)}$ and their corresponding parameters $\boldsymbol{\rm a}_{k}$ are generated by a beta process random measure and the conjugate Bernoulli process (BeP),

with base measure $B_{0}$ over the parameter space $\Theta=\mathbb{R}^{r}$ for our $r$-order autoregressive parameters $\boldsymbol{\rm a}_{k}$. As specified in Eq. (4), we take the normalized measure $B_{0}/B_{0}(\Theta)$ to be ${\cal N}(\boldsymbol{\rm m}_{0},\Sigma_{0})$. The discrete measures $B$ and $X^{(i)}$ can be represented as

with $f^{(i)}_{k}\sim\mbox{Ber}(\omega_{k})$.

The resulting feature vectors $\boldsymbol{\rm f}^{(i)}$ constrain the set of available states $z_{t}^{(i)}$ can take by constraining each transition distributions, $\boldsymbol{\rm\pi}^{(i)}_{j}$, to be 0 when $f_{k}^{(i)}=0$. Specifically, the BP-AR-HMM defines $\boldsymbol{\rm\pi}^{(i)}_{j}$ by introducing a set of gamma random variables, $\eta_{jk}^{(i)}$, and setting

The positive elements of $\boldsymbol{\rm\pi}^{(i)}_{j}$ can also be thought of as a sample from a finite Dirichlet distribution with only $K^{(i)}$ dimensions, where $K^{(i)}=\sum_{k}f^{(i)}_{k}$ represents the number of states channel $i$ uses. For convenience, we sometimes denote the set of transition variables $\{\eta_{jk}^{(i)}\}_{j}$ as $\boldsymbol{\rm\eta}^{(i)}$. As in the sticky HDP-HMM of Fox et al. [18], the parameter $\kappa_{c}$ encourages self-transitions (i.e., state $j$ at time $t-1$ to state $j$ at time $t$).

We again take a Bayesian nonparametric approach to define the event state HMM, building on the sticky HDP-HMM [18]. In particular, the transition distributions $\boldsymbol{\rm\phi}_{l}$ are hierarchically defined as

where $\mbox{stick}(\alpha)$ refers to a stick-breaking measure, also known as $\mbox{GEM}(\alpha)$, with $\boldsymbol{\rm\beta}$ generated by

Again, the sticky parameter $\kappa_{\rm e}$ promotes self-transitions, reducing state redundancy.

We term this model the sparse factorial BP-AR-HMM. Although the graph $G$ can be arbitrarily structured, because of our motivating seizure modeling application with a focus on a spatial-based graph structure, we often describe the sparse factorial BP-AR-HMM as capturing spatial correlations. We depict this model in the directed acyclic graphs shown in Figure 3. Note that while we formally consider a model of only a single event for notational simplicity, our formulation scales straightforwardly to multiple independent events. In this case, everything except the library of AR states $\{\boldsymbol{\rm a}_{k}\}$ becomes event-specific. If all events share the same channel setup, we can assume the channel covariances $\{\Delta_{l}\}$ are shared as well.

While interpreting these state sequences and covariances from the model, it is important to keep in mind that they are ultimately estimates of a system whose parsing even highly-trained physicians disagree upon. Nevertheless, we believe that the event states directly describe the activity of particular clinical interest.

In modeling the correlations between channels, the event states give insight into how different physiologic areas of the brain interact over the course of a seizure. In the clinical workup for resective brain surgery, these event states could help define and specifically quantify the full range of ways in which neurophysiologic regions initiate seizures and how others are recruited over the numerous seizures of a patient. In addition, given fixed model parameters, our model can fit the channel and event state sequences of an hour’s worth of 64-channel EEG data in a matter of minutes on a single 8-core machine, possibly facilitating epileptologist EEG annotation of long-term monitoring records.

The ultimate clinical aim of this work, however, involves understanding the relationship between epileptic bursts and seizures. Because the event state aspect of our model involves a Markov assumption, the intrinsic length of the event has little bearing on the states assigned to particular time points. Thus, these event states allow us to straightforwardly compare the neurophysiologic relationship dynamics in short bursts (often less than two seconds long) to those in much longer seizures (on the order of two minutes long), as explored in Section 4.4. Prior to this analysis, we first examine the importance of our various model components by comparing to baseline alternatives.

We explored the extent to which the spatial information and sparse dependencies encoded in the HIW prior improves our predictions of heldout data relative to a number of baseline models. To assess the impact of the sparse dependencies induced by the Gaussian graphical model for $\boldsymbol{\rm\epsilon}_{t}$, we compare to a full-covariance model with an IW prior on $\Delta_{l}$ (dense factorial). For assessing the importance of spatial correlations, we additionally compare to two alternatives where channels evolve independently: the BP-AR-HMM of Fox et al. [5] and an AR-HMM without the feature-based modeling provided by the beta process [24]. Both of these models use inverse gamma (IG) priors on the individual channel innovation variances. We learned a set of AR coefficients $\{\boldsymbol{\rm a}_{k}\}$ and event covariances $\{\Delta_{l}\}$ on one seizure and then computed the heldout log-likelihood on a separate seizure, constraining it use the learned model parameters from the training seizure.

For the training seizure, MCMC samples were collected over 5000 samples across 10 chains, each with a 1000-sample burn in and 10-sample thinning. To compute the predictive log-likelihood of the heldout seizure, we analytically marginalized the heldout event state sequence $Z_{1:T}$ but perform a Monte Carlo integration over the feature vectors $\boldsymbol{\rm f}^{(i)}$ and channel states $\boldsymbol{\rm z}_{1:T}$ using our MCMC sampler. For each original MCMC sample generated from the training seizure, a secondary chain is run fixing $\{\boldsymbol{\rm a}_{k}\}$ and $\{\Delta_{l}\}$ and sampling $z_{t}^{(i)}$, $Z_{t}$, $\boldsymbol{\rm f}^{(i)}$, $\boldsymbol{\rm\eta}^{(i)}$, and $\{\boldsymbol{\rm\phi}_{l}\}$ for the heldout seizure. We approximate $p(\boldsymbol{\rm y}_{1:T}\mid\{\boldsymbol{\rm\phi}_{l}\},\{\boldsymbol{\rm a}_{k}\},\{\Delta_{l}\})$ by averaging the secondary chain’s closed-form $p(\boldsymbol{\rm y}_{1:T}\mid\boldsymbol{\rm z}_{1:T},\{\boldsymbol{\rm\phi}_{l}\},\{\boldsymbol{\rm a}_{k}\},\{\Delta_{l}\})$, described in B.

Figure 6 (left) shows how conditioning on the innovations of neighboring channels in the sparse factorial model improves the prediction of an individual channel, as seen by its reduced innovation trace relative to the original BP-AR-HMM. The quantitative benefits of accounting for these correlations are seen in our predictions of heldout events, as depicted in Figure 6 (right), which compares the heldout log-likelihoods for the original and the factorial models listed above. As expected, the factorial models have significantly larger predictive power than the original models. Though hard to see due to the large factorial/original difference, the BP-based model also improves on the standard non-feature-based AR-HMM. Performance of the sparse factorial model is also at least as competitive as a full-covariance model (dense factorial). We would expect to see even larger gains for electrode grids with more channels due to the parsimonious representation presented by the graphical model. Regardless, these results demonstrate that the assumptions of sparsity in the channel dependencies do not adversely affect our performance.

In addition to providing a parsimonious modeling tool, the sparse dependencies among channels induced by the HIW prior allow our computations to scale linearly to the large number of channels present in iEEG. We compared a dense factorial BP-AR-HMM (entailing a fully-connected spatial graph) and a sparse factorial BP-AR-HMM on five datasets of 8, 16, 32, 64, and 96 channels (from three 32-channel electrodes) from the same seizure used previously. We ran the two models on each of the five datasets for at least 1000 MCMC iterations, using a profiler to tabulate the time spent in each step of the MCMC iteration.

Figure 7 shows the average time required to calculate the channel likelihoods at each time point under each AR channel state. This computation is used both for calculating the marginal likelihood (averaging over all the state sequences $z_{1:T}^{(i)}$) required in active feature sampling as well as in sampling the state sequences $z_{1:T}^{(i)}$. In our sparse factorial model, each channel has a constant set of $M$ dependencies, assuming $M$ neighboring channels. As such, the channel likelihood computation at each time point has an $O(M)$ complexity, implying an $O(MN)$ complexity for calculating the likelihoods of all $N$ channels at each time point. In contrast, the likelihood computation at each time point under the full covariance model had complexity $O(N)$, implying $O(N^{2})$ for calculating all the channel likelihoods. For $M\ll N$, as is typically the case, our sparse dependency model is significantly more computationally efficient.

Anecdotally, we also found that the IW prior experiments—especially those with larger number of channels—tended to occasionally have numerical underflow problems associated with the inverse term $\Delta^{-1(\boldsymbol{\rm i}^{\prime},\boldsymbol{\rm i}^{\prime})}_{Z_{t}}$ in the conditional channel likelihood calculation. This underflow in the IW prior model calculations is not surprising since the matrices inverted are of dimension $N-1$ (for $N$ channels), whereas in the HIW prior, the sparse spatial dependencies of the electrode grids make these matrices no larger than eight-by-eight.