Noisy network attractor models for transitions between EEG microstates

TL;DR

This paper introduces a novel stochastic differential equation model with noisy network attractors to simulate EEG microstate transitions, capturing their complex temporal dynamics and long-range correlations.

Contribution

It proposes a new nonlinear stochastic modeling framework that accurately reproduces EEG microstate transition probabilities and dwell-time distributions, including heavy tails.

Findings

The two-layer network model captures heavy-tailed dwell times.

Model fits EEG microstate sequences from different recording conditions.

Processing EEG data from MRI scanners affects dwell-time distribution tails.

Abstract

The brain is intrinsically organized into large-scale networks that constantly re-organize on multiple timescales, even when the brain is at rest. The timing of these dynamics is crucial for sensation, perception, cognition and ultimately consciousness, but the underlying dynamics governing the constant reorganization and switching between networks are not yet well understood. Functional magnetic resonance imaging (fMRI) and electroencephalography (EEG) provide anatomical and temporal information about the resting-state networks (RSNs), respectively. EEG microstates are brief periods of stable scalp topography, and four distinct configurations with characteristic switching patterns between them are reliably identified at rest. Microstates have been identified as the electrophysiological correlate of fMRI-defined RSNs, this link could be established because EEG microstate sequences are scale-free and have long-range temporal correlations. This property is crucial for any approach to model EEG microstates. This paper proposes a novel modeling approach for microstates: we consider nonlinear stochastic differential equations (SDEs) that exhibit a noisy network attractor between nodes that represent the microstates. Using a single layer network between four states, we can reproduce the transition probabilities between microstates but not the heavy tailed residence time distributions. Introducing a two layer network with a hidden layer gives the flexibility to capture these heavy tails and their long-range temporal correlations. We fit these models to capture the statistical properties of microstate sequences from EEG data recorded inside and outside the MRI scanner and show that the processing required to separate the EEG signal from the fMRI machine noise results in a loss of information which is reflected in differences in the long tail of the dwell-time distributions.

Figures

Tables

Table 1. Table 1: Two phase decay captures residence time distributions.

Table 2. Table 2: Parameters for the one-layer model.

Table 3. Table 3: Parameters for the two-layer model.

Equations

$$\mathcal{G}(t)=\sqrt{\frac{{\sum }_{i=1}^n(v_i(t)-\bar{v}(t){)}^2}{n}}$$

$$\sum _{i=1}^{k-1}\rho (i)<j\le \sum _{i=1}^k\rho (i)$$

$$T(m,j)=\frac{\#\{k:\sigma (k)=m\text{mbox}and\sigma (k+1)=j\}}{\#\{k:\sigma (k)=m\}}.$$

$$\begin{split}\tau\mathrm{d}p_{j}&=\left[f(p_{j},y)\right]+\eta_{p}\mathrm{d}w_{j}\\ \tau\mathrm{d}y_{k}&=\left[g\big{(}y_{k},A-Bp_{\alpha(k)}^{2}+C(y^{2}-y_{k}^{2})\big{)}\right]\mathrm{d}t+\eta_{{y}_{k}}\mathrm{d}w_{k}\end{split}$$

$$f(p_j,y)$$

$$g(y_k,\lambda )$$

$$p^2=\sum _{j=1}^M{p}_j^2,p^4=\sum _{j=1}^M{p}_j^4,y^2=\sum _{j=1}^Q{y}_j^2$$

$$\begin{array}{rl}{Z}_j^{(O)}(p,y)& =\sum _{\{k:\alpha (k)=j\}}{y}_k^2{p}_{\omega (k)}{p}_j\\ Z_j^{(I)}(p,y)& =\sum _{\{k^{\prime }:\omega (k^{\prime })=j\}}{y}_{k^{\prime }}^2{p}_{\alpha (k^{\prime })}^2.\end{array}$$

$$A=0.5,B=1.49,C=2,D=10,E=4,F=2$$

$$\begin{array}{rl}\tau \mathrm{d}{p}_{l,j}& =\left[f_l(p_j,y)\right]\mathrm{d}t+{\eta }_p\mathrm{d}w\\ \tau \mathrm{d}{y}_{l,k}& =\left[g_l(y_k,p_{\alpha (k)})+z_{l,k}\right]\mathrm{d}t+{\eta }_{l,k}\mathrm{d}{w}_k.\end{array}$$

$$z_{l+1,k}=\sum _{j=1}^{M_l}{\zeta }_{l,j}{p}_{l,j}^2.$$

$$f_l(p_j,y)$$

$$g_l(y_k,p_{\alpha (k)})$$

$$\begin{array}{rl}{Z}_{l,j}^{(O)}(p,y)& =\sum _{\{k:\alpha (k)=j\}}{y}_{l,k}^2{p}_{l,\omega (k)}{p}_{l,j}\\ Z_{l,j}^{(I)}(p,y)& =\sum _{\{k^{\prime }:\omega (k^{\prime })=j\}}{y}_{l,k^{\prime }}^2{p}_{l,\alpha (k^{\prime })}^2,\end{array}$$

$$p_l^2=\sum _{j=1}^{M_l}{p}_{l,j}^2,p_l^4=\sum _{j=1}^{M_l}{p}_{l,j}^4,y_l^2=\sum _{k=1}^{Q_l}{y}_{l,k}^2.$$

$${\mathcal{E}}^{(1)}(t)=A_1\exp (-k_{1}t),{\mathcal{E}}^{(2)}(t)=A_1\exp (-k_{1}t)+A_2\exp (-k_{2}t),$$

Full text

Introduction

The human brain is intrinsically organized into large-scale networks that can be identified at rest [1, 2]. These networks have to reorganize on a sub-second temporal scale in order to allow the precise execution of mental processes [3]. Spatial and temporal aspects of the dynamics underlying the reorganization of these large-scale networks require non-invasive measures with high spatial (functional magnetic resonance imaging (fMRI)) and temporal resolution (electroencephalography (EEG)). While fMRI uses the blood oxygenation level dependent (BOLD) response as a proxy for neuronal activity with a temporal resolution of several seconds, the EEG is a direct measure of neuronal activity which captures the temporal evolution of the scalp electrical field with millisecond resolution. Unlike local measures of the EEG in channel space that vary from time-point to time-point and as a function of the reference, the global measure of EEG topography remains stable for brief periods (50–100 ms) before changing to another quasi-stable state, the so-called EEG microstates [4, 5]. Interestingly, four dominant topographies are consistently reported both in healthy individuals as well as in neurological and psychiatric patients at rest [6, 7, 5]. While neurological and psychiatric diseases rarely affect their topography, they fundamentally alter their temporal dynamics [8, 9, 10]. They have been coined the “atoms of thought and can be considered the basic building blocks of spontaneous mentation that make up the spontaneous electrophysiological activity measured at the scalp [11].

A study using simultaneously recorded EEG-fMRI identified EEG microstates as the electrophysiological correlate of four fMRI-defined resting-state networks (RSNs)[12]. This link is surprising because EEG microstates and fMRI RSNs are two global measures of large-scale brain activity that are observed at temporal scales two orders of magnitude apart: 50–100 ms (microstates) and 10–20 seconds (fMRI-RSNs). Convolving the rapidly changing microstate time series with the hemodynamic response function (HRF) to obtain regressors for BOLD estimation acts like a strong temporal low-pass filter. The microstate time series is ’scale-free’ i.e. it shows the same behavior at different temporal scales, therefore no information is using this approach. Later, it was confirmed that this link could be established because EEG microstate time-series are mono-fractal and show long-range dependency (LRD) over six dyadic scales spanning two orders of magnitude (256 ms to 16 s) [13]. EEG microstates and fMRI RSNs hence capture the same underlying physiological process at different temporal scales with an electrophysiological and a hemodynamic measure, respectively. Importantly, Britz et al. demonstrate that the precise timing but not the order of local transitions of the microstate sequences is crucial for their fractal properties: shuffling their local transitions without changing their timing has no effect, whereas equalizing their durations degrades the time series to white noise without memory, hence attempts at modeling microstate sequences have to go beyond modelling the local transitions.

We note that the “state transition process” and the waiting (residence) time distribution for the “renewal process” where the transitions occur are essentially independent processes. In the case that the renewal process is memoryless then the whole process can be seen as Markov jump process, but it is possible for the transitions to be Markov but the jump times to be non-Markov. von Wegner et al. show that neither memoryless Markov models nor single parameter LRD models fully capture the data and conclude that more sophisticated models need to be developed to understand the underlying mechanisms of microstates [14].

In this paper we provide a novel modelling approach based on dynamical structures called noisy network attractors. These are stochastic models that exhibit heteroclinic or excitable network attractors in their noise-free dynamics [15]. A heteroclinic network is a collection of solutions (heteroclinic orbits) that link a set of steady states (saddles) that themselves are unstable. Excitable networks, in the sense introduced in [15], are close relations of heteroclinic networks that have a small but finite threshold of perturbation that needs to be overcome to make a transition between a number of attracting states. Heteroclinic networks have been found in models of many natural systems, for example from neuroscience [16, 17], population dynamics [18] and game theory [19]. The dynamics near a network attractor looks like a long transient: trajectories spend long periods of time close to one state before switching to another. Such transient dynamics have been observed for neural processes at a variety of levels of description [20] and examples ranging from olfactory processing in the zebrafish [21] to human working memory [22] have been successfully modeled using heteroclinic cycles or networks. Similar networks with noise have previously been shown to produce non-Markovian dynamics [23, 15]. Given the powerful capacities of network attractors to model transient dynamics at different levels these are a promising candidate to model EEG microstate sequences.

The transition probabilities and residence times in the EEG microstates sequence are modeled using SDEs that possess such a noisy network attractor, which is excitable (as described above) in the absence of noise. The noisy network attractor captures statistical properties of the observed microstate sequences including the distributions of residence times in each microstate and the transition probabilities. We apply this model to the analysis of EEG microstate sequences obtained from resting state EEG recordings reported in [13]. We show that the distributions of times spent in each microstate best fit a sum of exponentials, and so the transitions are non Markov. We construct one and two layer models and apply each model to microstate sequences. We demonstrate that the one layer model produces a single exponential distribution of residence times with the same local transition probabilities as the data but no LRD. We further show that the double layer model is required to achieve distributions that are a sum of exponentials that capture both the local transition probabilities, distribution of dwell times and LRD and provides hence a better fit to the data.

Methods

Data collection

Detailed description of the procedures used to collect the EEG recordings and convert them into microstates are given in [13]; here we provide a brief summary for completeness.

Subjects and Procedure

Nine healthy volunteers (24 – 33 years, mean age 28.37 years) participated for monetary compensation after giving informed consent approved by the ethics commission of the University Hospital of Geneva. None suffered from current or prior neurological or psychiatric illness or from claustrophobia. For each subject, we recorded one 5-minute session outside the scanner prior to recording three 5-minute runs of resting-state EEG inside the MRI scanner. Subjects were instructed to relax and rest with their eyes closed without falling asleep and to move as little as possible. Data from one subject had to be excluded due to subsequent self-report of sleep and the presence of sleep patterns in the EEG, and the data from the remaining eight subjects were submitted to further analysis.

EEG Recording

The EEG was recorded from 64 sintered Ag/AgCL electrodes mounted in an elastic cap and arranged in an extended 10-10 system. The electrodes were equipped with an additional 5$k\Omega$ resistor, and impedances were kept below 15$k\Omega$. The EEG was digitized using a battery-powered and MRI-compatible EEG system (BrainAmp MR plus, Brainproducts) with a sampling frequency of 5 kHz and a hardware bandpass filter of 0.016 – 250 Hz with the midline fronto-central electrode as the physical reference. The amplifier was placed ca. 15 cm outside the magnet bore and data were transmitted via fiberoptic cables to the recording computer placed outside the scanner room.

EEG Data Preprocessing

Three types of artifacts were removed for data recorded inside the scanner: first, gradient artifacts were removed using a sliding average [24] and then the EEG was down-sampled to 500 Hz and low-pass filtered with a finite-impulse response filter with a low-pass of 70 Hz. Next, the ballistocardiographic (BCG) artifact was removed using a sliding average and finally, independent component analysis (ICA) was used to remove the residual BCG artifact along with oculomotor and myogenic artifacts. Data recorded outside the scanner were first downsampled to 500 Hz and subsequently, oculo-motor and myogenic artifacts were removed using ICA. Finally, both the EEG recorded inside and outside the scanner was further downsampled to 125 Hz and bandpass filtered between 1 and 40 Hz.

EEG Microstates

The Global Field Power (GFP) is a measure of the overall strength of the scalp electrical field. Between the local troughs of the GFP, the scalp topography remains stable and only varies in strength. Hence, the local peaks of the GFP are the best representative of an EEG microstate. We extracted the EEG at all local peaks of the GFP and submitted those to a modified Atomize-Agglomerate Hierarchical (AAHC) clustering method [25] in order to determine the best solution representing the most dominant topographies. The best solution was identified using a cross-validation criterion [26], a measure of residual variance, which identified four template maps as the optimal solution in each run. Finally, we computed the spatial correlation between the four dominant template maps of the optimal solution and the continuous EEG data using a temporal constraint criterion of 32 ms to obtain the time-course of the dominant EEG microstates.

In total we process and analyze 8 recordings outside the scanner and 24 recordings inside the scanner. We convert the resting state EEG into time series where each time point of the recording is classified into one of the four classically identified microstates [6, 7, 5]. The raw data is first cleaned to remove oculo-motor and other artifacts, and resampled at 125Hz. For each time point of the cleaned data we computed the global field potential

[TABLE]

where $v_{i}(t)$ is the voltage of channel $i$ and $\bar{v}(t)$ is the voltage mean. The EEG topography at the local maxima of the $\mathcal{G}$ has a good signal to noise ratio and so we take these points as momentary maps. All momentary maps then undergo a clustering analysis to identify the four dominant map topographies, these are the microstates. The four microstates are shown in Figure 1B. Finally, each time point of the (resampled) recording is classified into one of the four possible microstates.

EEG Microstate Sequence Analysis

We denote the sequence of microstates $m(t)$ where $m\in\{1,2,3,4\}$ at any given sampling time point $t$. We classify $m(t)$ into epochs of the same state defined by saying that we enter a new epoch if and only if $t=0$ or $m(t+1)\neq m(t)$.

Analogous to the model output, we describe the $k$th epoch in terms of its state $\sigma(k)$ and residence time $\rho(k)$ respectively. This means we choose the unique sequence such that if

[TABLE]

then $m(j)=\sigma(k)$ and moreover $\sigma(k)\neq\sigma(k+1)$ (there is a change of state at the end of each epoch).

Specifically, we define

• •

$R(t)$ is the distribution of residence times $\rho(k)$ for all epochs $k$.

• •

$T(m,j)$ is the probability of transition from an epoch in state $m$ to one in state $j$,

[TABLE]

We call the sequence of nodes visited $\sigma(k)$ the transition process, and a sequence of residence times $\rho(k)$ the renewal process. These processes are essentially independent. The state transition process is often represented as a transition matrix containing the probabilities $T$. Each probability only depends on the state that you are currently in. In the case that the renewal process is memoryless, the whole process can be seen as Markov jump process, but it is possible for the transitions to be Markov but the renewal process to be non-Markov. In the latter case the distribution of residence times would be not be a simple exponential.

The average number of epochs per recording for outside is 2010, and for inside is 3031. We compute $R$ by plotting a histogram of the residence times for each recording, then find the average and standard error of each bin over all recordings in that group (outside or inside). Similarly we calculate $T$ for each recording and find the average and standard error for each probability.

Mathematical models

Single-layer network model

We aim to build a model that captures the statistical properties of the transition and renewal processes. To this end, we construct a model of stochastic differential equations perturbed by low amplitude additive white noise, using a general method that allows us to realize any desired graph as an attracting heteroclinic or excitable network in phase space. This model has evolved from work in [27] and is detailed in [15] so we briefly outline the construction.

We realize the network of all possible transitions between the four canonical microstates as an excitable network with four nodes ($p$-cells). There is an edge between two nodes in the network if there can be transition between the two corresponding microstates; here the network is all-to-all connected with twelve edges ($y$-cells). The coupling architecture of the network and corresponding microstates are shown in Figure 1A–B.

The system is given by:

[TABLE]

for $j=1,\cdots,M$ and $k=1,\cdots,Q$. The functions $f$ and $g$ are defined by

[TABLE]

where

[TABLE]

and the outputs $(O)$ and inputs $(I)$ to the $p$ cells from the $y$ cells are:

[TABLE]

The $w$ are independent identically distributed noise processes, $\eta$ are noise weights and $A,B,C,D,E,F$ are constants. Here we introduce the time scaling $\tau$ because although the $p$ cells can be scaled by the parameters the $y$ cells have a functional form which is fixed. The full equations are given in S1 Appendix..

The $p_{j}$ variables classify which node of the network (i.e. which of the four microstates) is visited. In the system with no noise, there are equilibrium solutions with $p_{j}=1$ and all other coordinates zero. The $y_{k}$ variables become non-zero during a transition between nodes.

We use the parameter set

[TABLE]

throughout as in [15], where $B<1.5$ gives an excitable network. In an excitable network nodes are stable equilibria and transitions between nodes are driven by the additive noise. The transition rates between nodes are modulated by the noise levels on the edges (rather than on the nodes) as described in [15]. We choose the noise levels so that the model statistics, namely the transition probabilities the distribution of the renewal process (residence times) fit the data; see Model simulation and fitting for details. We fix the noise on the nodes $\eta_{p}=10^{-4}$. The time scaling constant $\tau$ is set to 8ms (in line with the sampling rate of the data) throughout.

Figures 1C–D show example time series output for the nodes and the edges, respectively, for the excitable network model. At any given moment in the time series the trajectory in the simulation is close to one of the equilibrium solutions, and one of the $p_{j}$ variables is close to 1. To determine where the trajectory is close to a given node we define a box in phase space around that node so that when the trajectory is in the box we say it is near the node. Here we fix the box size $H=0.49$ such that each box contains one node and the boxes do not overlap. The duration of time spent in the box is then the residence time for that node [15].

The output of the model simulation is a series of $k$ epochs where each epoch is defined by the state it is in $\sigma\in\{1,2,3,4\}$ and its residence time $\rho$. The transition and renewal processes generated by the network construction given here are both Markov. Due to the evidence of multi-scale behavior of the temporal dynamics of microstate sequences [28, 13] we also present a more sophisticated model a with a multi-layer construction that will generate a non-Markov renewal process.

Multi-layer network model

We construct a system of $N$ levels, where each level $l$ has $M_{l}$ nodes, and we assume all-to-all connections, so we thus have $Q_{l}\equiv M_{l}(M_{l}-1)$ edges. In each level, we set up a system of SDEs in the form of that described in [15], as follows, where $p_{l,j}\in\mathbb{R}^{M_{l}}$ and $y_{l,k}\in\mathbb{R}^{Q_{l}}$:

[TABLE]

Here we allow for a general input into the $y$-cells $z_{l,k}(t)$ that linearly couples layers from the $p_{l}$ nodes to the $y_{l+1,k}$ edges by:

[TABLE]

The parameter $\zeta_{l,j}$ is a constant that scales the activity from node $p_{l,j}$ to all the edges in level $l+1$. The functions $f$, $g$, $Z^{(O)}$ and $Z^{(I)}$ are analogous to those in the single layer model,

[TABLE]

for $j=1,\dots,M_{l}$ and $k=1,\dots,Q_{l}$,

[TABLE]

and

[TABLE]

We consider in particular a two layer model where $M_{1}=2$ and $M_{2}=4$ and the constants $A,B,C,D,E,F$ are set as (5) for each layer. Note that both layers are excitable networks. Layer one with two nodes is a “hidden layer" and we only consider the output from layer 2 with four nodes. As before the output is a transition process and a renewal process.

Figure 2A shows the coupling architecture for the two layer model, B–C show the time series of the nodes and edges in layer one, respectively, and D–E show the time series of the nodes and edges in layer two, respectively. Compare layer two topography and dynamics to Figure 1. The two nodes in layer one effect the dynamics on the edges (and therefore the residence times) in layer two by the scalings $\zeta_{1,1}=\zeta_{1}$ and $\zeta_{1,2}=\zeta_{2}$. These scale the dwell times in the renewal process output from layer two. For illustrative purposes we choose $\zeta_{1}=10^{-1}$ and $\zeta_{2}=10^{-3}$ here. Panels B–E clearly show that when $p_{1,2}=1$ the residence times at each node $p_{2,j}$ are longer (transitions between nodes are less frequent) as the edges in layer two $y_{2,k}$ are scaled by $\zeta_{2}$. Note that if $\zeta_{1}=\zeta_{2}$ layer one would be redundant as the residence times would be consistent (drawn from the same distribution) in either node and the renewal process would again be Markov.

Model simulation and fitting

We wish to use the model capture the transition probabilities and distribution of residence times of the data. To this end we first fit exponential curves to the residence distributions for the data recorded inside and outside the fMRI scanner. Specifically, we compute the histogram of residence times from the microstate sequences for bin size $40$. The distribution is truncated at the first empty bin. We then use the fit function in Matlab to fit the single and double exponentials

[TABLE]

to each dataset where we constrain $A_{i}>0$ and $k_{i}\geq 0$. We use an F-test to indicate whether the double exponential was a better fit to the data; code used as written for [29].

We numerically simulate the one and two-layer models with a stochastic Heun method implemented using a custom code written in Matlab. We compute 10 realizations of the model using step size $0.05$ up to maximum of 100,000 steps. This gives approximately 5750 epochs per realization. We calculate the residence time distribution $R$ and transition probabilities $T$, then average the values over 10 realizations. The results of the averaged simulations (model output) are compared by eye to the data.

To fit the one-layer model output to the data we adjust the noise amplitude parameters on the edges $\eta_{y_{k}}$. The noise amplitude on each edge controls the probability of transition along that edge, therefore we change the amplitudes so the model output transition probabilities are within the error bars of the transition probabilities from the data. For example if $T(1,2)>T(2,1)$ then we set $\eta_{y_{1}}>\eta_{y_{4}}$. The overall magnitude of the $\eta_{y_{k}}$ parameters controls the decay rate of the residence time distribution. If all $\eta_{y_{k}}$ are large, for example, $O(10^{-1})$ transitions happen quickly, the residence times are short, and the decay rate of the distribution $R$ is large (slope is steep); whereas if $\eta_{y_{k}}$ are small, for example, $O(10^{-4})$ there are long residence times between transitions and the decay rate is small (slope is shallow). In this way we can fit the distribution of the model to the fitted curves ${\mathcal{E}}^{(1)}(t)$.

To fit the transition probabilities in the two-layer model we adjust the noise values of the edges in layer two $\eta_{2,k}$ in the same way as for the one-layer model. To fit the residence distributions we first fix the overall magnitude of the $\eta_{2,k}$ in line with the one-layer model, then adjust the transfer parameters $\zeta_{j}$ and noise values of the edges on layer one $\eta_{1,k}$. Changing $\zeta_{j}$ changes the dynamics on the edges of layer two in a homogeneous way (the same scaling on all edges): if $\zeta$ is increased the residence times decrease and the decay rate of the distribution is becomes larger; if $\zeta$ is decreased, the decay rate of the distribution of the residence times is decreased. As there are two transfer parameters (one for each node in layer one) the residence distribution of the output from the model is a linear combination of the two decay rates; one $\zeta$ is associated with the distribution at short times and the other captures the heavy tails. The proportion of each distribution (the mixing) is controlled by the noise on the edges is layer one. If $\eta_{1,1}>\eta_{1,2}$ more time will be spent in node $p_{1,2}$ than in node $p_{1,1}$ so in the residence distribution there will be a larger proportion of the decay rate associated with $\zeta_{2}$. In this way we can fit the distribution of the model to the fitted curves ${\mathcal{E}}^{(2)}(t)$.

Results

Microstate residence time distributions have two phase decay

The residence time distributions $R(t)$ for data collected inside and outside the scanner and the transition probabilities $T(m,j)$ for the EEG recordings are shown in Figure 3. The distributions $R(t)$ are plotted on a logarithmic scale with ${\mathcal{E}}^{(1)}(t)$ and ${\mathcal{E}}^{(2)}(t)$ given by (8). The transition probabilities for each dataset are shown in panels C.

For each dataset we perform a comparison of fits of ${\mathcal{E}}^{(1)}$ and ${\mathcal{E}}^{(2)}$ to the data. We find that an $F$-test rejects the null hypothesis that the distribution is ${\mathcal{E}}^{(1)}$: details are given in Table 1. Using the Kruskall-Wallis ANOVA the mean rank fo the residence time distributions are significantly different ($p<0.0005$). Significant differences in transition probabilities between outside and inside can be seen between in the transition from state 1 to state 2, and from state 4 to state 2. Note there are no self transitions due to our definition of an epoch.

One-layer network model captures transition probabilities

The results of fitting the one-layer model to the data from outside the scanner are shown in Figure 4. The single exponential fit ${\mathcal{E}}^{(1)}$ from Table 1 is shown. The transition probabilities for the model and the data show a good fit as the model probabilities are within the error bars of the data for all transitions. The noise weights $\eta_{y_{k}}$ used are given in Table 2. The distribution of the one-layer model aligns with ${\mathcal{E}}^{(1)}$. However, this model fails to capture the two decay rates observed in the data.

Two-layer model captures two-phase decay of residence times

The results of fitting the two-layer model to the data from outside the scanner are shown in Figure 5 with the double exponential fit ${\mathcal{E}}^{(2)}$ from Table 1. The distributions of the simulations agree closely with ${\mathcal{E}}^{(2)}$ and fit within the error bars of the distribution. The transition probabilities for the model and data are also shown. The parameters for the simulations are shown in Table 3. The two layer model captures the two-phase decay of the data.

Discussion

In this article we have demonstrated a new modelling approach applied to capture the statistical properties of the temporal dynamics of EEG microstates. We analyze EEG microstate sequences from EEG data recorded inside and outside an fMRI scanner, previously reported in [12, 13]. We consider the transition probabilities between microstates and the distributions of the residence times.

We show that there are significant differences between the residence time distributions. Both residence time distributions are best fit by two decay rates. We note the similarity in the first decay rates given by $2.340\times 10^{-2}$ and $2.620\times 10^{-2}$ for outside and inside respectively. However, the second decay rate is very different $6.249\times 10^{-3}$ and $1.198\times 10^{-2}$ for outside and inside respectively. We also see a change in the transition probabilities; the inside dataset transition probabilities have less variance than the outside dataset, with a significant increase of the probability of transition to microstate 2 from microstates 1 and 4. The most likely explanation for the difference in long-range temporal correlation for the data recorded inside and outside the scanner is the difference in artifact removal. The data recorded outside the scanner were only contaminated by very occasional oculomotor artifacts (eye-movements) which can be reliably detected and removed with ICA. EEG data recorded inside the scanner are additionally contaminated by the gradient and the ballistocardiographic artifacts. The gradient artifact has an amplitude three orders of magnitude larger than the EEG, but it is very regular and can thus be removed reliably by a sliding average. The amplitude of the ballistocardiogram is roughly an order of magnitude larger than the EEG, and its removal requires two steps: a sliding average to remove the major amplitude and ICA to remove the residual. The latter is essential but also prone to remove not only the residual artifacts but also to remove part of the signal.

In the past few years there have been various attempts to model EEG-microstate sequences. Gärtner et al. [30] construct a hidden-Markov type stochastic model based on the transition probabilities between microstates extracted from the data. The transition matrix gives the probabilities of moving from one state to another, where each probability only depends on the state that you are currently in. Such a one-step Markov construction assumes that the transitions between microstates depend on the current state but otherwise are memoryless. However, no restrictions are placed on the distributions of the residence times.

The Markov model approach taken in [30] has been criticized for the underlying assumption that the microstate transition process is independent of the underlying global field power time series [31]. Further, the authors do not comment on possible long range dependencies (LRDs) in the data. Gschwind et al. argue that LRDs are an intrinsic property of neural dynamics that manifest in the microstate sequences as scale-free or fractal structure in the data [28]. To elucidate this fractal structure and verify and expand the work of [13] they subject the data to a battery of tests, including computing the power spectral density and Hurst exponents using detrended fluctuation analysis (DFA), a wavelet framework, and time-variance analysis. Using wavelet detrended fluctuation analysis EEG-microstate time series have been shown to exhibit mono-fractal behavior over six dyadic scales indicative of long range  [12].

Further analysis of scale-free properties in EEG-microstate data, subsequently conducted by von Wegner et al., shows that these measures, when applied to real (instead of simulated) EEG-data, do not guarantee the existence of LRDs [14]. Specifically, they compute the Hurst exponent of a given EEG-microstates data set using three different methods and perform statistical tests that show that the results are not always significantly different from a constructed white noise process. The assumptions underlying the analysis that provide a connection between the Hurst exponent and long range dependencies in data and therefore temporal fractal structure, may also not be applicable to real EEG-data.

That the microstate residence times presented here are best represented by two decays rates provides additional evidence that the temporal properties of microstates are non-Markovian in line with [28]. It is clear that Markovian or single parameter LRD models do not adequately represent the temporal dynamics of microstate sequences [14].

One limitation of this approach is that we only consider memoryless transitions, and assume that the transition probabilities do not change over time. These assumption may be incorrect as they do not fit with the non-Markovian properties of microstate sequences previously found [13]. Future work would be to look at transition probabilities of sequences of microstates, thereby including different levels of history.

Combined fMRI-EEG studies have found correlates between each microstate and resting state networks [13]; A is the auditory network, B is the visual network, C is the saliency network and D is the attention network. The switching dynamics between microstates indicates changes in switching between these underlying neural networks [12]. Multiple decay rates point to complex switching dynamics required for cognition. More work needs to be done to understand how these dynamic changes correlate to the cognitive dysfunction observed in this state, for example loss of short term memory.

Finally this modelling approach could be used to identify dynamic differences and their brain network correlates at different levels of consciousness and in other neurological disorders for example Alzheimer’s disease, epilepsy and schizophrenia.

Conclusions

Using EEG data from [13] we show that the residence time distributions of resting state EEG microstate sequences have a two phase decay. We also show that the additional processing of the data recorded inside the fMRI scanner leads to a decrease of very long residence times ($>900$ms) in the distribution. Using the construction outlined in [15] we build systems of stochastic differential equations (SDE) that have a “noisy network attractor”. Using a one-layer model we show that this captures the transition probabilities between microstates. The distribution given by this model is (single) exponential. We also build a two-layer model with a (hidden) layer containing two nodes each associated with a transfer parameter that scales the behavior of the edges in the four-node network layer. This model produces a residence time distribution with two decay rates (sum of two exponentials) that better captures the two-phase decay seen in the data. The identification of EEG-microstates as correlates of resting state brain networks means that analysis of the temporal dynamics of microstate sequences can lead to crucial insights into the switching between these networks required for cognition and perception in resting state and could provide indicators of underlying mechanisms of dynamic changes characteristic of neurological disorders.

Supporting information

S1 Appendix.

Single-layer model equations.

Acknowledgements

JC and PA gratefully acknowledge the financial support of the EPSRC Centre for Predictive Modelling in Healthcare, via grant EP/N014391/1. JC acknowledges funding from MRC Skills Development Fellowship MR/S019499/1. CMP acknowledges funding from the Marsden Fund, Royal Society of New Zealand.