IDENTIFICATION OF DYNAMIC FUNCTIONAL BRAIN NETWORK STATES THROUGH
TENSOR DECOMPOSITION

arXiv:1410.0446v1 [] 2 Oct 2014

Arash Golibagh Mahyari, Selin Aviyente
mahyari@msu.edu, aviyente@egr.msu.edu
Department of Electrical and Computer Engineering
Michigan State University
East Lansing, MI 48823 USA
ABSTRACT
With the advances in high resolution neuroimaging, there has
been a growing interest in the detection of functional brain
connectivity. Complex network theory has been proposed as
an attractive mathematical representation of functional brain
networks. However, most of the current studies of functional
brain networks have focused on the computation of graph theoretic indices for static networks, i.e. long-time averages of
connectivity networks. It is well-known that functional connectivity is a dynamic process and the construction and reorganization of the networks is key to understanding human
cognition. Therefore, there is a growing need to track dynamic functional brain networks and identify time intervals
over which the network is quasi-stationary. In this paper, we
present a tensor decomposition based method to identify temporally invariant ’network states’ and find a common topographic representation for each state. The proposed methods
are applied to electroencephalogram (EEG) data during the
study of error-related negativity (ERN).
Index Terms— Graphs, Dynamic Networks, Tensor Decomposition, Electroencephalography.
1. INTRODUCTION
With the advance in noninvasive imaging modalities such
as fMRI, EEG, and MEG, it is important to develop computational methods capable of giving a succinct description
of the functional brain networks [1]. Functional connectivity describes the coordinated activation of segregated groups
of neurons. Traditionally, functional connectivity has been
quantified through linear [2] and nonlinear measures [3].
Synchronization of neuronal oscillations has been suggested
as one plausible mechanism in the interaction of spatially
distributed neural populations and has been quantified using
phase synchrony [4]. Although phase synchrony is successful
at quantifying pairwise interactions [5], it cannot completely
This work was in part supported by the National Science Foundation
under Grant No. CCF-1218377.

describe the complex relationship between function and organization of the brain. Recently, research in the area of
complex networks has led to fundamental insights into the
organization of the healthy and diseased brain [5, 6]. However, the current studies are limited to the analysis of static
brain networks obtained through averaging long-term functional connectivity and thus, neglect possible time-varying
properties of the topologies.
There is growing evidence that functional networks dynamically reorganize and coordinate on millisecond scale for
the execution of mental processes [7]. For this reason, there
has been an interest in characterizing the dynamics of functional networks using high temporal resolution EEG recordings. The early work in this area was an extension of static
network analysis to the dynamic case by extracting graph theoretic features from graphs across time and tracking the evolution of these parameters [8]. However, these approaches
lose the spatial information provided by the graphs and cannot identify which parts of the brain contributed to the observed changes in the network. More recently, a ”network
state” framework has been proposed [9, 7], where each state is
defined as periods of time during which the network topology
is quasi-stationary. In this paper, we adopt this framework for
tracking the topology of brain networks across time and representing each state with a common topographic map. The current work differs from the existing approaches in a couple of
ways. First, we take into account the full structure of the network at each time point. Second, we consider extracting network states common across time and subjects unlike current
work which considers individual subjects. Finally, the current work offers a compressed spatial representation of each
network state through tensor-tensor projection unlike current
approaches which use averaging. Tensor-to-tensor projection
proposed in this paper projects the information across subjects
and time into a lower dimensional ’signal’ subspace whereas
averaging assigns equal weights to all subjects.
It is also important to note that the proposed framework is
closely tied to dynamic network tracking. The most common
approaches to network tracking have been to identify anomalies using subspace projection methods such as in [10, 11] or

through sliding window estimates with time independent or
dependent weighting factors [12, 13]. More recently, adaptive
evolutionary clustering [14] was proposed to track the cluster
changes over time. However, all of these methods either detect time points where events of interest happen, or find different clustering structures at each time point. In this paper, we
propose a comprehensive framework that first identifies time
intervals during which the network topology is stationary and
then summarizes each interval by a single lower dimensional
network through tensor-tensor projection.
2. BACKGROUND
2.1. Time-Varying Network Construction
The time-varying functional brain networks are constructed
from multichannel EEG data with the nodes corresponding
to the different brain regions and the edges to the connectivity between these regions. In this paper, we quantify the
connectivity using a recently introduced phase synchrony
measure based on RID-Rihaczek distribution [15]. The first
step in quantifying phase synchrony is to estimate the time
and hfrequency
i dependent phase, Φi (t, ω), of a signal, si ,
arg

Ci (t,ω)
|Ci (t,ω)|
1

Therefore, {G(t)}t=1,2,...,T is a time series of N × N
weighted and undirected graphs corresponding to the functional connectivity network at time t for a fixed frequency
band, where T is the total number of time points and N is the
number of nodes within the network.
2.2. Tensor Subspace Analysis
Linear data models such as Principal Component Analysis (PCA) and Independent Component Analysis (ICA) are
widely used for the decomposition of matrices. Depending
on the criteria, different types of basis vectors are extracted
and appropriate lower dimensional features are determined
through projection. Multiway data analysis extends these linear methods to capture multilinear structures and underlying
correlations in higher-order datasets, also known as tensors.
Some exemplary methods include PARAFAC, Tucker decomposition, and Higher-Order Singular Value Decomposition
(HOSVD) [16, 17].
The Tucker decomposition is a higher order generalization of Singular Value Decomposition (SVD) [16]. Let X ∈
Rm1 ×m2 ...×md be a d -mode array, then its Tucker decomposition can be expressed as:

, where Ci (t, ω) is the complex RID-Rihaczek

distribution :
Ci (t, ω) =

Z Z


θτ
(θτ )2
exp(j )Ai (θ, τ )e−j(θt+τ ω) dτ dθ
exp −
σ
2
|
{z
} | {z }


Choi-Williams kernel

Rihaczek kernel

(1)

R
and Ai (θ, τ ) = si (u + τ2 )s∗i (u − τ2 )ejθu du is the ambiguity
function of the signal s i . The phase synchrony between nodes
i and j at time t and frequency ω is computed using Phase
Locking Value (PLV):
P LVij (t, ω) =

L

1 X
exp jΦkij (t, ω)
L

where L is the number of trials and Φkij (t, ω) = |Φki (t, ω) −
Φkj (t, ω)| is the phase difference between the two channels for
the k th trial.
Once the pairwise synchrony values are computed at each
time and frequency point, we can construct a time-varying
graph {G(t)}t=1,2,...,T with
ωb
1 X
P LVij (t, ω)
Ω ω=ω

X = C ×1 U(1)×2 U(2) ... ×d U(d)+ E
(1)

(2)

(d)

Ci1 ,i2 ,...,id ui1 ◦ ui2 ... ◦ uid

+ Ei1 ,i2 ,...,id
(4)
where C ∈ Rr1 ×r2 ...×rd is the core tensor, and U(1) ∈
Rm1 ×r1 , U(2) ∈ Rm2 ×r2 , ..., U(d) ∈ Rmd ×rd , where
r1 ≤ m1 , r2 ≤ m2 , ..., rd ≤ md , are the projection matrices whose columns are orthogonal. E ∈ Rm1 ×m2 ...×md
is the residual error, and ×k is the product of a tensor and a
matrix along mode-k. Reconstruction of the original tensor,
X̃ ∈ Rm1 ×m2 ...×md , from a lower dimensional representation is obtained as:
i1 ,i2 ,...,id

(2)

k=1

Gij (t) =

=

P

(3)

a

where Gij (t) ∈ [0, 1] represents the connectivity strength between the nodes i and j within the frequency band of interest,
[ωa , ωb ], and Ω is the number of frequency bins in that band.
1 The details of the RID-Rihaczek distribution and the corresponding synchrony measure are given in [15].







†
†
†
X˜ = X ×1 U(1) U(1) ×2 U(2) U(2) ... ×d U(d) U(d)
(5)
where † is the transpose of the matrix.
3. THE PROPOSED METHOD
3.1. Temporal Tracking for Network State Identification
In the proposed work, the time-varying functional connectivity graphs across subjects are considered as a 4-mode tensor
G ∈ RN ×N ×T ×S constructed as channel × channel × time
× subject, with N being the number of channels, T the total number of time points and S the number of subjects. The
Tucker decomposition of this connectivity tensor yields:
G = C ×1 U(1) ×2 U(2) ×3 U(3) ×4 U(4) + E

(6)

where C ∈ RN ×N ×T ×S is the core tensor, and E N ×N ×T ×S
is the residual error.
To obtain an approximation of G, G̃ ∈ RN ×N ×T ×S , we
first consider the full Tucker decomposition with U(1) ∈
RN ×N , U(2) ∈ RN ×N , U(3) ∈ RT ×T , and U(4) ∈ RS×S .
The singular values along each mode are ordered by fixing
the index of all of the other modes to 1. Since first singular
values along each mode represent the largest variance of the
data along that mode, we choose that to order the remaining mode. To get the approximation tensor, the appropriate
number of singular vectors along first and second modes N̄
is defined as N̄ = jk , where jk is the highest index for which
|Cjk ,1,1,1 | ≥ 0. Similarly, the number of singular vectors
along the fourth mode is S̄ = sk , where sk is the highest
index for which |C1,1,1,sk | ≥ 0. The time mode is not projected to a lower dimensional space since all time points are
necessary to identify the exact boundaries of the network
states. The lower dimensional projection matrices are de(2) (2)
(2)
(1) (1)
(1)
fined as: Ũ(1) = [u1 u2 ...uÑ ], Ũ(2) = [u1 u2 ...uÑ ],
(4) (4)

(4)

Ũ(4) = [u1 u2 ...uS̃ ].
The reconstructed tensor G̃ ∈ RN ×N ×T ×S is obtained as:




†
†
G̃ = G ×1 Ũ(1) Ũ(1) ×2 Ũ(2) Ũ(2)




†
†
×3 Ũ(3) Ũ(3) ×4 Ũ(4) Ũ(4) .

(7)

The 4-mode approximation tensor G̃ can be written as a
sequence of 3-mode tensors G̃t ; t = 1, 2, ..., T . To detect
the boundaries of network states, we propose a new temporal clustering algorithm. Unlike regular data clustering, the
proposed method considers both the similarity of the lower
dimensional representation of the networks as well as their
closeness in time. The similarity of two networks at time
t1 and t2 is quantified through a cosine similarity metric
between G̃t1 and G̃t2 as follows:
∆(t1 , t2 ) =

hG̃t1 ,G̃t2 i
kG̃t1 kkG̃t2 k

; t1 , t2 = 1, 2, ..., T

(8)

where ha, bi is the inner product of a and b, and k a k is the
Frobenius Norm. Similarly, the temporal closeness between
−(t1 −t2 )2

two graphs is quantified as Θ(t1 , t2 ) = e 2σ2 ; t1 , t2 =
1, 2, ..., T , where σ is a parameter which determines the
weighting for different time separations, and depends on the
sampling frequency.
The combined similarity matrix is defined as:
Ψ(t1 , t2 ) = λΘ(t1 , t2 ) + (1 − λ)∆(t1 , t2 ); t1 , t2 = 1, 2, ..., T
(9)
where λ ∈ (0, 1) determines the trade-off between tensor
similarity and time proximity. This similarity matrix is input to a standard spectral clustering algorithm combined with
k-means to identify the boundaries of the network states [18].

3.2. Topographic Compression for Network State Representation
Once the time boundaries of the different network states are
identified, each state has to be summarized with a single topographic map. Previously, this was commonly addressed by
averaging the edges over the time interval [19]. This method
has the drawback of emphasizing all of the edges equally and
resulting in very dense network representations.
For a given time interval (T1 , T2 ) and the 3-mode tensor
sequence corresponding to this interval GT1 , GT1 +1 , GT1 +2 , ..., GT2 ,
the goal is to extract the topographic map Ĝ ∈ RN ×N which
best represents that network state. The 3-mode tensors corresponding to the time interval (T1 , T2 ) can be rewritten
as a 4-mode tensor by taking the time modality into account G ′ ∈ RN ×N ×(T2 −T1 +1)×S , and decomposed using
(1)
(2)
the full Tucker decomposition, U′
∈ RN ×N , U′
∈
(3)
(4)
RN ×N , U′ ∈ R(T2 −T1 +1)×(T2 −T1 +1) , U′ ∈ RS×S , similar to Equation 6.
In order to summarize the subject information to find a
general unique model which fits all subjects, the 4-mode tensor G ′ ∈ RN ×N ×(T2 −T1 +1)×S is projected by the singular
(4)
vector u′l corresponding to the lth largest singular value in
this mode. Likewise, to summarize the time information of
the resulting 3-mode tensor, it is projected to the singular vec(3)
tor u′k corresponding to the k th largest singular value of the
(3)
(4)
time mode, Ĝ = G ′ ×3 u′k ×4 u′l . The values of k and l
are usually equal to 1 but may change depending on the data.
4. EXPERIMENTAL RESULTS
4.1. EEG Data
The proposed framework is applied to a set of EEG data
containing the error-related negativity (ERN) 2 . The ERN is
a brain potential response that occurs following performance
errors in a speeded reaction time task usually 25-75 ms after
the response [20]. Previous work [21] indicates that there
is increased coordination between the lateral prefrontal cortex (lPFC) and medial prefrontal cortex (mPFC) within the
theta frequency band (4-8 Hz) and ERN time window (25-75
ms), supporting the idea that frontal and central electrodes
are functionally integrated during error processing. EEG
data from 62-channels was collected in accordance with the
10/20 system on a Neuroscan Synamps2 system (Neuroscan,
Inc.). A speeded-response flanker task was employed, and
response-locked averages were computed for each subject.
All EEG epochs were converted to current source density
(CSD) using published methods [22]. Data were averaged
across trials (∼ 200 trials) for the purpose of ERN and timefrequency analysis.In this paper, we analyzed data from 91
subjects corresponding to the error responses.
2 We thank Dr. Edward Bernat from the University of Maryland for sharing his EEG dataset with us.

4.2. Network State Identification and Summarization
The connectivity matrices are constructed by computing the
pairwise average PLV between 62 channels in the theta frequency band for all time (2 seconds) and all subjects using
Eq. 2. The time-varying graphs {G(t)}t=1,2,...,T for all subjects and all time will be treated as a 4-mode tensor, which is
decomposed using Tucker decompostion. The approximation
tensors G̃t ; t = 1, 2, ..., T with N̄ = 2, S̄ = 3 are used to
obtain the 256 × 256 similarity matrix Ψ. The matrix Ψ is
computed with λ = 0.4 and σ = 2500 as shown in Fig. 1.
The values of λ and σ are empirically chosen to obtain the
best separation between clusters. Once the matrix Ψ is ob-

long-range relationships we are interested in.
In order to show the most significant edges in the summarized graph, the edges with values in the top 1% are selected
and plotted in Fig. 2. Most of the significant edges are in the
frontal and central areas where node AF 8 acts as a hub, with
the highest degree equal to 59. C1 has the second highest
degree (53) followed by F P 2 (52). These nodes correspond
to the right lateral prefrontal cortex (lPFC) and medial prefrontal cortex (mPFC) in accordance with previous findings
which indicate increased synchronization between these regions during ERN [21].

60

FP1 FPZ FP2
AF8
AF7
AF3 AFZ
F7 F5
F8
F1
F3
FZ F2 F4 F6
FT8
FT7FC5FC3
FC1FCZFC2FC4FC6
T7 C5 C3 C1 CZ C2 C4 C6 T8
CPZCP2CP4CP6
TP7CP5CP3CP1
TP8
PZ
P2 P4 P6
P1
P7 P5 P3
P8
PO3 POZPO4
PO6
PO5
PO8
PO7
O1 OZ O2

50

40

30

20

10

Fig. 2. The topographic map with the most signifcant edges
of the compressed connectivity matrix for the network state
(−132ms, 188ms).
Fig. 1. Similarity matrix, Ψ, computed for 2 seconds of EEG
data across 91 subjects in the theta frequency band.
5. CONCLUSIONS
tained, the critical time points are detected using the spectral
clustering with K = 5. The number of clusters K is selected
based on the eigenspectrum of the similarity matrix [23]. The
detected time intervals are (−1000, −703) ms, (−703, −132)
ms, (−132, 188) ms,(188, 736) ms,(736, 1000) ms.
As expected, the first two time intervals correspond to the
prestimulus part where there is less change in the network
configuration. The third time interval (−132ms, 188ms) is
of particular interest since it includes the time interval right
before a response is made as well as the ERN interval (0 −
150ms). This is extracted as a separate network state indicating a reorganization of the functional network configuration.
Similarly, the time interval (188, 736)ms contains the P300
event which is expected to result in a distinct topographic
map. In this paper, we will focus on extracting the common
topographic map for the time interval (−132, 188) ms since
it coincides with the ERN [21]. To obtain a single network
representation for this time interval, we selected k = l = 2
instead of the singular vectors corresponding to the highest
singular values since the projection to the subspace spanned
by the largest singular value mostly contains edges between
physically adjacent nodes. This is a side effect of volume conduction affecting PLV values and does not convey the actual

In this paper, we proposed a tensor-based method for monitoring dynamic functional connectivity networks to identify
quasi-stationary network states and represent the common topographic distribution of each state. Network states were
identified computing a similarity metric which takes both the
similarity of the reconstructed tensors and their proximity in
time to obtain a partitioning of the dynamic networks into
contiguous time intervals. After identifying the boundaries of
the network states, a topographical map for each time interval
was obtained by a tensor-tensor projection. The application
of the proposed algorithm to ERN data yields time intervals
that closely correspond to events of interest and topographic
maps that are consistent with previous hypotheses with regard
to error monitoring in the brain. The proposed method is time
consuming for large data sets due to the complexity of Tucker
decomposition.
Future work will consider extensions of this framework
to partitioning time and frequency dependent connectivity
networks by considering higher order tensor representations.
Moreover, the choice of optimal parameters (λ and σ) will
be considered using cost functions such as modularity for
evaluating the quality of the different partitions.

6. REFERENCES
[1] B. He, T. Coleman, G.M. Genin, G. Glover, X. Hu,
N. Johnson, T. Liu, S. Makeig, P. Sajda, and K. Ye,
“Grand challenges in mapping the human brain: NSF
workshop report,” IEEE transactions on Biomedical engineering, 2013. in press.
[2] M. Rubinov and O. Sporns, “Complex network measures of brain connectivity: uses and interpretations,”
Neuroimage, vol. 52, no. 3, pp. 1059–1069, 2010.
[3] E. Pereda, R.Q. Quiroga, and J. Bhattacharya, “Nonlinear multivariate analysis of neurophysiological signals,”
Progress in Neurobiology, vol. 77, no. 1-2, pp. 1–37,
2005.
[4] F. Varela, J.P. Lachaux, E. Rodriguez, and J. Martinerie,
“The brainweb: Phase synchronization and large-scale
integration,” Nature Reviews Neuroscience, vol. 2, no.
4, pp. 229–239, 2001.
[5] E. Bullmore and O. Sporns, “Complex brain networks:
graph theoretical analysis of structural and functional
systems,” Nature Reviews Neuroscience, vol. 10, no.
3, pp. 186–198, 2009.
[6] C.J. Stam and J.C. Reijneveld, “Graph theoretical analysis of complex networks in the brain,” Nonlinear
biomedical physics, vol. 1, no. 1, pp. 3, 2007.
[7] S.I. Dimitriadis, N.A. Laskaris, and A. Tzelepi, “On
the quantization of time-varying phase synchrony patterns into distinct functional connectivity microstates
(fcµstates) in a multi-trial visual ERP paradigm,” Brain
topography, pp. 1–13, 2013.
[8] S.I. Dimitriadis, N.A. Laskaris, V. Tsirka, M. Vourkas,
S. Micheloyannis, and S. Fotopoulos, “Tracking brain
dynamics via time-dependent network analysis,” Journal of Neuroscience Methods, vol. 193, no. 1, pp. 145–
155, 2010.
[9] R.F. Betzel, M.A. Erickson, M. Abell, B.F. O’Donnell,
W.P. Hetrick, and O. Sporns, “Synchronization dynamics and evidence for a repertoire of network states in
resting eeg,” Frontiers in computational neuroscience,
vol. 6, 2012.
[10] A. Lakhina, M. Crovella, and C. Diot, “Diagnosing
network-wide traffic anomalies,” in ACM SIGCOMM
Computer Communication Review. ACM, 2004, vol. 34,
pp. 219–230.
[11] B.A. Miller, M.S. Beard, and N.T. Bliss, “Eigenspace
analysis for threat detection in social networks,” in 2011
Proceedings of the 14th International Conference on Information Fusion (FUSION), 2011, pp. 1–7.

[12] H. Tong, S. Papadimitriou, P. S. Yu, and C. Faloutsos,
“Proximity tracking on time-evolving bipartite graphs,”
in 8th SIAM International Conference on Data Mining,
Applied Mathematics.
[13] K. Xu, M. Kliger, and A. Hero, “A shrinkage approach
to tracking dynamic networks,” in Proceedings of the
IEEE Statistical Signal Processing Workshop (SSP ’11),
2011, p. 517520.
[14] K.S. Xu, M. Kliger, and A.O. Hero III, “Adaptive evolutionary clustering,” arXiv preprint arXiv:1104.1990,
2011.
[15] S. Aviyente and A.Y. Mutlu, “A time-frequency-based
approach to phase and phase synchrony estimation,”
IEEE Transactions on Signal Processing, vol. 59, no.
7, pp. 3086–3098, 2011.
[16] E. Acar and B. Yener, “Unsupervised multiway data
analysis: A literature survey,” Knowledge and Data Engineering, IEEE Transactions on, vol. 21, no. 1, pp. 6–
20, 2009.
[17] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A
multilinear singular value decomposition,” SIAM journal on Matrix Analysis and Applications, vol. 21, no. 4,
pp. 1253–1278, 2000.
[18] U. Von Luxburg, “A tutorial on spectral clustering,”
Statistics and computing, vol. 17, no. 4, pp. 395–416,
2007.
[19] A. G. Mahyari and S. Aviyente, “Two-dimensional
svd for event detection in dynamic functional brain networks,” in Proceedings of IEEE Global Conference on
Signal and Information Processing, 2013.
[20] J. R. Hall, E. M. Bernat, and C. J. Patrick, “Externalizing psychopathology and the error-related negativity,” Psychological Science, vol. 18, no. 4, pp. 326–333,
2007.
[21] J.F. Cavanagh, M.X. Cohen, and J.J.B. Allen, “Prelude
to and resolution of an error: EEG phase synchrony reveals cognitive control dynamics during action monitoring,” The Journal of Neuroscience, vol. 29, no. 1, pp.
98–105, 2009.
[22] J. Kayser and C.E. Tenke, “Principal components analysis of laplacian waveforms as a generic method for identifying ERP generator patterns: I. evaluation with auditory oddball tasks,” Clinical Neurophysiology, vol. 117,
no. 2, pp. 348–368, 2006.
[23] A.Y. Ng, M.I. Jordan, and Y. Weiss, “On spectral clustering: Analysis and an algorithm,” Advances in neural information processing systems, vol. 2, pp. 849–856,
2002.

