Low-Dimensional Dynamics of Resting-State Cortical Activity

Abstract

Endogenous brain activity supports spontaneous human thought and shapes perception and behavior. Connectivity-based analyses of endogenous, or resting-state, functional magnetic resonance imaging (fMRI) data have revealed the existence of a small number of robust networks which have a rich spatial structure. Yet the temporal information within fMRI data is limited, motivating the complementary analysis of electrophysiological recordings such as electroencephalography (EEG). Here we provide a novel method based on multivariate time–frequency interdependence to reconstruct the principal resting-state network dynamics in human EEG data. The stability of network expression across subjects is assessed using resampling techniques. We report the presence of seven robust networks, with distinct topographic organizations and high frequency (∼5–45 Hz) fingerprints, nested within slow temporal sequences that build up and decay over several orders of magnitude. Interestingly, all seven networks are expressed concurrently during these slow dynamics, although there is a temporal asymmetry in the pattern of their formation and dissolution. These analyses uncover the complex temporal character of endogenous cortical fluctuations and, in particular, offer an opportunity to reconstruct the low dimensional linear subspace in which they unfold.

Appendices

Appendix 1: Multivariate Decomposition

Interdependences between channel pairs were assessed in the frequency domain. Spectral decomposition of EEG signals were obtained using the Fourier transform of sliding time windows as

$$ \begin{aligned} S_{xy}[n,k]&=S_x[n,k]S^{*}_y[n,k]=|S_{x}[n,k]||S_{y}[n,k]|\hbox{e}^{{\bf j}[\phi_x[n,k]-\phi_y[n,k ] ] },\\ S_{x}[n,k]&=\sum_{l=1}^L x_n[l]w_n[l]e^{-{\bf j}\frac{2\pi}{L}kl}, \end{aligned} $$

(4)

where ${\bf j}=\sqrt{-1}, \phi_x[n,k]$ and ϕ _y [n, k] denote the phase-spectra, S _x [n, k] denotes the short-time Fourier transform of a discrete-time signal x[p] in the nth time-window, w _n [l]. The w[l] is the unit energy window function and L is the total number of data samples in each time-window. The asterisk (*) denotes the complex conjugate. Time–frequency interdependence was assessed between all channel pairs. The total number of channel pairs is given by

$$ B=\frac{c(c-1)}{2}=1830, $$

(5)

where c = 61 denotes the number of EEG channels.

A single multidimensional array for each subject $j=1,2, \ldots ,7$ and epoch i = 1, 2, 3 was constructed as

$$ \hat{\Upxi}_{i,j_{(r,s)}}[n,k]= \left[\begin{array}{ccccc} \hat{\theta}_{I_{(1,2)}}[n,k] &\hat{\theta}_{I_{(1,3)}}[n,k] & \hat{\theta}_{I_{(1,4)}}[n,k] & \cdots & \hat{\theta}_{I_{(1,61)}}[n,k] \\ \dagger & \hat{\theta}_{I_{(2,3)}}[n,k] & \hat{\theta}_{I_{(2,4)}}[n,k] & \cdots & \hat{\theta}_{I_{(2,61)}}[n,k]\\ \dagger & \dagger & \hat{\theta}_{I_{(3,4)}}[n,k] & \cdots & \hat{\theta}_{I_{(3,61)}}[n,k]\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \dagger & \dagger & \cdots & \dagger & \hat{\theta}_{I_{(60,61)}}[n,k] \end{array}\right], $$

(6)

where $\hat{\theta}_{I_{(r,s)}}[n,k]$ denotes the imaginary part of the time–frequency interdependence expressed in Eq. (1) between channels r and s. The symbol $\dagger$ denotes array elements that, due to the undirected nature of Eq. (1) are symmetric across the major diagonal and these are not considered further. Since Eq. (1) is a lossless decomposition of the EEG data, the elements in the array $\hat{\Upxi}_{i,j_{(r,s)}}[n,k]$ will inevitably be mutually correlated and contain considerable redundancies.

PCA was implemented through an eigen-decomposition of the covariance matrix of the array $\hat{\Upxi}_{i,j_{(r,s)}}[n,k].$ Since $\hat{\Upxi}_{i,j_{(r,s)}}[n,k]$ has three factors (or dimensions) of time (n), frequency (k), and channel-combinations (b), we reshaped $\hat{\Upxi}_{i,j_{(r,s)}}[n,k]$ into a 2-D array $\Uptheta_{i,j}[nk,b], \quad nk=1,2, \ldots ,N\times K; \quad b=1,2, \ldots ,B,$ where B denotes the maximum number of possible channel combinations as defined in Eq. (5). That is, we reorganized the dimensions of $\hat{\Upxi}_{i,j_{(r,s)}}[n,k]$ as time–frequency (N × K)-by-channel-combinations (B), expressed in matrix format $\Uptheta_{i,j_{(nk\times n)}}.$ We then obtained the resting-state modes from the eigen-decomposition of the covariance matrix of $\Uptheta_{i,j}[nk,b]$ as

$$ \hat{R}_{\Uptheta\Uptheta_{i,j}}={{\mathbb{E}}}\{\Uptheta_{i,j}^T\Uptheta_{i,j}\} =E_{i,j}^{T}D_{i,j}, $$

(7)

where ${\mathbb{E}}$ denotes the mathematical expectation, T the matrix transpose and, as above, i denotes each 4-min epoch within each subject j. E _i,j is the B × B orthogonal matrix whose columns are the unit-norm eigenvectors of the covariance matrix $\hat{R}_{\Uptheta\Uptheta_{i,j}},$ and D _i,j the diagonal matrix representing the eigenvalues of $\hat{R}_{\Uptheta\Uptheta_{i,j}},$ such that $\hbox{diag}\Big(D_{i,j}\Big)=\{d_{i,j}^1>d_{i,j}^2> \ldots >d_{i,j}^B\}$ represents the variance of each putative resting-state mode. The eigenvectors of $\hat{R}_{\Uptheta\Uptheta_{i,j}}$ correspond to the edges of the resting-state networks for the ith epoch in the jth subject. That is, each B × 1 eigenvector contains the contributions of every channel combination (i.e., edges) to that specific mode and hence reflects the spatial topology of each network. We considered the first 100 modes (columns of E _i,j ) for further analysis. After realignment, this yields 100 modes for each epoch i and subject j as

$$ \begin{array}{l} Z_{i,j}[nk,m]=\Uptheta_{i,j}[nk,b]E_{i,j}[b,m]\\ m=1,2, \ldots ,B; \quad m=1,2, \ldots ,100 . \end{array} $$

(8)

Jackknifing revealed that seven resting-state networks were robustly expressed across subjects and epochs. These networks were considered for further analysis. Although the variance between epochs can be used to assess the inter-subject (or test-retest) reliability, we focus on between-subject reliability and hence average across epochs, which yields

$$ \begin{aligned} {\bf E}[b,m]&=\left[E_{j}[b,m]\right],\\ {\bf D}[m]&=\left[D_{j}[m,m]\right],\\ {\bf Z}[nk,m]&=\left[Z_{j}[nk,m]\right], \end{aligned} $$

(9)

where, as before $j=1,2, \ldots ,7$ runs across the seven subjects and $m=1,2, \ldots ,7$ denotes the seven robust resting-state networks. Hence E[b, m]_{(B × 7 × 7)}, D[m]_(7 × 7), and Z[nk,m]_{(NK × 7 × 7)} represent the 7 robust resting-state modes (eigenvectors), explained variances (eigenvalues), and projection time–frequency spectra for 7 subjects, respectively. For example, $Z_{1}[nk,2]=\Uptheta_{1}[nk,b]E_{1}[b,2]$ is the projection time–frequency spectra for the second network in subject 1. The matrix E[b, m] contains the functional connectivity structure of each mode that can be visualized in terms of nodes (i.e., EEG electrode sites) or edges (i.e., the weighted links between nodes). After reshaping, the time–frequency spectra Z _j,m [n, k] contain the spectral content and the temporal fluctuations of each robust networks.

Appendix 2: Temporal Expression of Robust Networks

We further assessed the networks’ temporal dynamics by estimating the Hurst exponent $\hat{H}_{j,m}=\frac{\hbox{log}\it{(R_{j,m}/S_{j,m})}}{\hbox{log}(N)},$ where R _j,m /S _j,m denotes the so-called range statistics: R _j,m is the difference between the maximum and minimum deviation from the mean, and S _j,m is the standard deviation over total number of N samples for the jth subject and the mth network. The Hurst exponent is a measure of the extent of long-range dependence, i.e., for $\hat{H}_{j,m}=0.5$ increments that are uncorrelated, whereas $0.5<\hat{H}_{j,m}<1$ indicates the presence of long-range dependencies or persistence (Kobayashi et al. 1999).

The Hurst exponent of the temporal fluctuations $\acute{\sigma}^2_{j,m}[n]$ of each network expressed in Eq. (3), is given by

$$ \frac{R_{j,m}}{S_{j,m}}=\frac{\underset{1\leq M\leq N}{\max}\vartheta_{j,m}-\underset{1\leq M\leq N}{\min}\vartheta_{j,m}} {\sqrt{\frac{1}{N}\sum_{n=1}^{N}(\acute{\sigma}^2_{j,m}[n] -\bar{\acute{\sigma}}_{j,m})^2}}, $$

(10)

where

$$ \begin{aligned} \vartheta_{j,m}&=\sum_{l=1}^{M} \acute{\sigma}^{2}_{j,m}[l]-\bar{\acute{\sigma}}_{j,m},\\ \bar{\acute{\sigma}}_{j,m}&=\frac{1}{N} \sum_{n=1}^{N}\acute{\sigma}^2_{j,m}[n], \end{aligned} $$

(11)

where ϑ _j,m denotes the deviation of mth network from its mean value $\bar{\acute{\sigma}}_{j,m},$ and non-overlapped samples are chosen as M = N/4.

We finally study the relative temporal expression of resting-state networks by determining the cross-correlation between their time series $\acute{\sigma}^{2}_{j,m}[n].$ To ensure normality, the time-series $\acute{\sigma}^{2}_{j,m}[n]$ were first log transformed and then low-pass filtered at frequency 0.1 Hz. The cross-correlation between the times series of two networks of u and v is given as a function of a variable time lag τ by

$$ \begin{aligned} C_{j_{(u,v)}}[\tau]&={{\mathbb{E}}}\{\acute{\sigma}^{2}_{j,u}[n] \acute{\sigma}^{2}_{j,v}[n+\tau]\}, \quad u\neq v=1,2, \ldots ,7,\\ \tau&=-N,-N+1, \ldots ,0,1,2, \ldots ,N. \end{aligned} $$

(12)

To extract the component of the cross-correlation functions that is not invariant to time reversal, we determine the asymmetric components of C _{j_(u,v)}[τ] by reconstructing

$$ \acute{C}_{j_{(u,v)}}[\tau]=\frac{C_{j_{(u,v)}}[\tau] -C^{T}_{j_{(u,v)}}[\tau]}{2}, $$

(13)

where C _{j_(u,v)}[τ] is the cross-correlation function of the time series $\acute{\sigma}^{2}_{j,m}[n]$ and T the transpose. The signs of the slope of the asymmetric component at zero lag, $\acute{C}_{j_{(u,v)}}[0],$ were used to represent the relative leading (positive sign) and lagging (negative sign) position of networks u and v. Examining all pairs of networks in each subject allows a temporal rank order (sequence) to be obtained.

Confidence intervals for the asymmetric component of the cross correlation were estimated using a non-parametric permutation approach. To this end, 1000 surrogate time series were constructed by random permutation of the log transformed low-pass filtered network time series. The asymmetric component of the cross-correlation functions Eq. (13) derived from these surrogate time series hence represent trivial non-zero fluctuations of the magnitude expected from data of that string length and amplitude distribution but with no further structure of interest. Two-sided 95 % confidence intervals were estimated by rank ordering this null distribution and choosing the 25th and 975th values. This process was repeated independently for each of the seven robust networks in all 3 epochs of all subjects.