1

Graph Frequency Analysis of Brain Signals

arXiv:1512.00037v2 [q-bio.NC] 3 May 2016

Weiyu Huang, Leah Goldsberry, Nicholas F. Wymbs, Scott T. Grafton, Danielle S. Bassett and Alejandro Ribeiro

Abstract—This paper presents methods to analyze functional
brain networks and signals from graph spectral perspectives.
The notion of frequency and filters traditionally defined for
signals supported on regular domains such as discrete time and
image grids has been recently generalized to irregular graph
domains, and defines brain graph frequencies associated with
different levels of spatial smoothness across the brain regions.
Brain network frequency also enables the decomposition of brain
signals into pieces corresponding to smooth or rapid variations.
We relate graph frequency with principal component analysis
when the networks of interest denote functional connectivity. The
methods are utilized to analyze brain networks and signals as
subjects master a simple motor skill. We observe that brain
signals corresponding to different graph frequencies exhibit
different levels of adaptability throughout learning. Further, we
notice a strong association between graph spectral properties of
brain networks and the level of exposure to tasks performed,
and recognize the most contributing and important frequency
signatures at different task familiarity.
Index Terms—Functional brain network, network theory,
graph signal processing, fMRI, motor learning, filtering

I. I NTRODUCTION
The study of brain activity patterns has proven valuable
in identifying neurological disease and individual behavioral
traits [1]–[4]. The use of functional brain networks describing
the tendency of different regions to act in unison has proven
complementary in the analysis of similar matters [5]–[8]. It is
not surprising that signals and networks prove useful in similar
problems since the two are closely related. In this paper we
advocate an intermediate path in which we interpret brain activity as a signal supported on the graph of brain connectivity.
We show how the use of graph signal processing (GSP) tools
can be used to glean information from brain signals using the
network as an aid to identify patterns of interest. The benefits
of incorporating network information into signal analysis has
been demonstrated in multiple domains. Notable examples of
applications include video compression [9], rating predictions
in recommendation systems [10], and breast cancer diagnostics
[11], and semi-supervised learning [12].
The fundamental GSP concepts that we utilize to exploit
brain connectivity in the analysis of brain signals are the
graph Fourier transform (GFT) and the corresponding notions
Supported by NSF CCF-1217963, PHS NS44393, ARO ICB W911NF09-0001, ARL W911NF-10-2-0022, ARO W911NF-14-1-0679, NIH R01HD086888, NSF BCS-1441502, NSF BCS-1430087 and the John D. and
Catherine T. MacArthur and Alfred P. Sloan foundations. Content does not
necessarily represent official views of any of the funding agencies. W. Huang,
L. Goldsberry, D. S. Bassett, and A. Ribeiro are with the Dept. of Electrical
and Systems Eng., University of Pennsylvania. N. F. Wymbs is with the
Dept. of Physical Medicine and Rehabilitation, Johns Hopkins University.
S. T. Grafton is with the Dept. of Psychological and Brain Sciences,
University of California at Santa Barbara. Email: whuang@seas.upenn.edu,
goleah@seas.upenn.edu, nwymbs1@jhmi.edu, scott.grafton@psych.ucsb.edu,
dsb@seas.upenn.edu, and aribeiro@seas.upenn.edu.

of graph frequency components and graph filters. These concepts are generalizations of the Fourier transform, frequency
components, and filters that are used in regular domains such
as time and spatial grids [13]–[15]. As such, they permit
the decomposition of a graph signal into components that
represent different modes of variability. We can define low
graph frequency components representing signals that change
slowly with respect to brain connectivity networks in a well
defined sense and high graph frequency components representing signals that change fast in the same sense. This is
important because low and high temporal variability have
proven important in the analysis of neurological disease and
behavior [16], [17]. GFT based decompositions permit a
similar analysis of variability across regions of the brain for a
fixed time – a sort of spatial variability measured with respect
to the connectivity pattern. We demonstrate here that it is
useful in a similar sense; see e.g. Figs. 5, 6, 7, and 9.
The GSP studies in this paper are related to principal
component analysis (PCA), which has been used with success
in the analysis of brain signals [18], [19]. The difference
with the GSP analysis we present here is that PCA implicitly
assumes the brain network to be a correlation matrix and the
signals to be drawn from a stochastic model. More importantly,
whereas the GFT can be used to, e.g., decompose the signal
into low, medium, and high frequency components, PCA is
mostly utilized for dimensionality reduction; which in the
language of this paper is tantamount to analyzing a few low
graph frequency components. Another important difference is
that PCA focuses on identifying variability across different
realizations of brain signals, but the GFT identifies spatial
variability of a single realization. GSP is also related to the
spectral analysis of networks in general and Laplacians in
particular [20], [21]. The difference in this case is that these
spectral analyses yield properties of the networks. In GSP
analyses, the network provides an underlying structure, but
the interest is on signals expressed on this stratum.
The goal of this paper is to introduce GSP notions that can
be used to analyze brain signals and to demonstrate their value
in identifying patterns that appear when monitoring activity as
subjects learn to perform a visual-motor task. Specifically, the
contributions of this paper are: (i) To explain tools from the
emerging field of GSP and show how they can be applied in
analyzing brain signals. (ii) To evaluate the graph spectrum
of brain functional network and to define artificial network
construction methods that replicate the features of the graph
spectrum of functional networks. (iii) To examine the temporal
variation of brain signals corresponding to different graph
frequencies when participants perform visual-motor learning
tasks. (iv) To investigate the contribution of brain signals
associated with different graph frequencies to the learning
success at different stages of visual-motor learning.

2

We begin the paper with the introduction of basic notions
of graphs and graph signals. Particular emphasis goes into the
definition of the graph Fourier transform and the interpretation
of graph frequency components as different modes of spatial
variability measured with respect to the brain network (Section
II-A). We also introduce the notion of graph filters and
discuss the interpretation of a graph low-pass filter as a local
averaging operation (Section II-B). Band-pass and high-pass
filters that are used in later sections to separate signals into
different components are also introduced. We point out that the
discussion here is more extensive than necessary for readers
familiar with GSP. The goal is to make the paper accessible
to readers that are not necessarily familiar with the subject.
We then move on to describe two different experiments
involving the learning of different visual-motor tasks by different sets of participants (Section III). The graph frequency
of functional brain networks for the participants is visualized
and analyzed (Section IV), which is found to be associated
across scan sessions in the same dataset and across datasets.
In specific, we find that graph high frequencies of functional
networks concentrate on visual and sensorimotor modules of
the brain – the two brain areas well-known to be associated
with motor learning [22], [23]. This motivates us to consider graph frequencies besides low frequency components,
whereas the PCA-oriented approaches has been focusing on
low frequencies. We also describe the construction of a simple
model to establish artificial networks with a few network
descriptive parameters (Section IV-A). We observe that the
model is able to mimic the properties of actual functional brain
networks and we use them to analyze spectral properties of the
brain networks (Section IV-B). The paper then utilizes graph
frequency decomposition to visualize and investigate brain
activities with different levels of spatial variation (Section V).
It is noticed that the decomposed signals associated to different
graph frequencies exhibit different levels of temporal variation
throughout learning (Section V-A). In particular, graph low and
high frequency components exhibit higher temporal variation
at multiple temporal scales in the two experiments considered.
Because temporal variation has been shown to be associated
with learning success [17], this implies graph low and high frequency components offer more contributions to learning success. Finally, we also define learning capabilities of subjects,
and examine the importance of brain frequencies at different
task familiarity by evaluating their respective correlation with
learning performance at different task familiarities (Section
VI). We find that it is good to have graph low frequency
components (smooth, spread, and cooperative brain signals)
when we face an unfamiliar task. When we become highly
familiar with the task, it is better and favors further learning
to have graph high frequency components (varied, spiking, and
competitive brain signals).
II. G RAPH S IGNAL P ROCESSING
The interest of this paper is to study brain signals in which
we are given a collection of measurements xi associated with
each cortical region out of n different brain regions. An
example signal of this type is an fMRI reading in which xi

estimates the level of activity of brain region i. The collection
of n measurements is henceforth grouped in the vector signal
x = [x1 , x2 , . . . , xn ]T ∈ Rn . A fundamental feature of
the signal x is the existence of an underlying pattern of
structural or functional connectivity that couples the values of
the signal x at different brain regions. Irrespective of whether
connectivity is functional or structural, our goal here is to
describe tools that utilize this underlying brain network to
analyze patterns in the neurophysiological signal x.
We do so by modeling connectivity between brain regions
with a network that is connected, weighted, and symmetric.
Formally, we define a network as the pair G = (V, W), where
V = {1, . . . , n} is a set of n vertices or nodes representing
individual brain regions and W ∈ Rn×n collects weights of
edges in the network with wij ≥ 0 the weight of the edge
(i, j). Since the network is undirected and symmetric we have
that wij = wji for all (i, j). The weights wij = wji represent
the strength of the connection between regions i and j, or,
equivalently, the proximity or similarity between nodes i and
j. In terms of the signal x, this means that when the weight
wij is large, the signal values xi and xj tend to be related.
Conversely, when the weight wij is small, the signal values
xi and xj are not directly related except for what is implied
by their separate connections to other nodes.
We adopt the conventional definitions of the degree and
Laplacian matrices [24, Chapter 1]. The degree matrix D ∈
n×n
R
P+ is a diagonal matrix with its ith diagonal element Dii =
j wij denoting the sum of all the weights of edges out of
node i. The Laplacian matrix is defined as the difference L :=
D − W ∈ Rn×n . The components of the Laplacian
matrix are
Pn
explicitly given by Lij = −wij and Lii = j=1 wij . Observe
that the Laplacian is real, symmetric, diagonal dominant, and
with strictly positive diagonal elements. As such, the matrix L
is positive semidefinite. The eigenvector decomposition of L
is utilized in the following section to define the graph Fourier
transform and the associated notion of graph frequencies.
We note that brain networks, irrespective of whether their
connectivity is functional [25] or structural [26], tend to be
stable for a window of time, entailing associations between
brain regions during captured time of interest. Brain activities
can vary more frequently, forming multiple samples of brain
signals supported on a common underlying network.
A. Graph Fourier Transform and Graph Frequencies
Begin by considering the set {λk }k=0,1,...,n−1 of eigenvalues of the Laplacian L and assume they are ordered so that
0 = λ0 ≤ λ1 ≤ . . . ≤ λn−1 . Define the diagonal eigenvalue
matrix Λ := diag(λ0 , λ1 . . . , λn−1 ) and the eigenvector matrix
V := [v0 , v1 , . . . , vn−1 ].

(1)

Because the graph Laplacian L is real symmetric, it accepts
the eigenvalue decomposition
L = VΛVH ,

(2)

where VH represents the Hermitian (conjugate transpose)
of the matrix V. The validity of (2) follows because the
eigenvectors of symmetric matrices are orthogonal so that the

3

definition in (1) implies that VH V = I. The eigenvector
matrix V is used to define the Graph Fourier Transform of
the graph signal x as we formally state next; see, e.g., [15].
Definition 1 Given a signal x ∈ Rn and a graph Laplacian
L ∈ Rn×n accepting the decomposition in (2), the Graph
Fourier Transform (GFT) of x with respect to L is the signal
x̃ = VH x.

(3)

The inverse (i)GFT of x̃ with respect to L is defined as
x = Vx̃.

(4)

We say that x and x̃ form a GFT transform pair.
Observe that since VVH = I, the iGFT is, indeed, the
inverse of the GFT. Given a signal x we can compute the
GFT as per (3). Given the transform x̃ we can recover the
original signal x through the iGFT transform in (4).
There are several reasons that justify the association of the
GFT with the Fourier transform. Mathematically, it is just a
matter of definition that if the vectors vk in (2) are of the form
vk = [1, ej2πk/n , . . . , ej2πk(n−1)/n ]T , the GFT and iGFT in
Definition 1 reduce to the conventional time domain Fourier
and inverse Fourier transforms. More deeply, it is not difficult
to see that if the graph G is a cycle, the vectors vk in (2)
are of the form vk = [1, ej2πk/n , . . . , ej2πk(n−1)/n ]T . Since
cycle graphs are representations of discrete periodic signals,
it follows that the GFT of a time signal is equivalent to the
conventional discrete Fourier transform; see, e.g., [27].
An important property of the GFT is that it encodes a notion
of variability akin to the notion of variability that the Fourier
transform encodes for temporal signals. To see this, define
x̃ = [x̃0 , . . . , x̃n−1 ]T and expand the matrix product in (4) to
express the original signal x as
x=

n−1
X

x̃k vk .

(5)

k=0

It follows from (5) that the iGFT allows us to write the
signal x as a sum of orthogonal components vk in which
the contribution of vk to the signal x is the GFT component x̃k . In conventional Fourier analysis, the eigenvectors
vk = [1, ej2πk/n , . . . , ej2πk(n−1)/n ]T carry a specific notion
of variability encoded in the notion of frequency. When k is
close to zero, the corresponding complex exponential eigenvectors are smooth. When k is close to n, the eigenfunctions
fluctuate more rapidly in the discrete temporal domain. In
the graph setting, the graph Laplacian eigenvectors provide a
similar notion of frequency. Indeed, define the total variability
of the graph signal x with respect to the Laplacian L as
X
TV(x) = xH Lx =
wij (xi − xj )2 ,
(6)
i6=j

where in the second equality we expanded the quadratic form.
It follows that the total variation TV(x) is a measure of how
much the signal changes with respect to the network. For the
edge (i, j), when wij is large we expect the values xi and xj
to be similar because a large weight wij is encoding functional
similarity between brain regions i and j. The contribution of

their difference (xi − xj )2 to the total variation is amplified
by the weight wij . If the weight wij is small, activities at
brain regions i and j tend to be uncorrelated, and therefore
the difference between the signal values xi and xj makes little
contribution to the total variation. We can then think of a signal
with small total variation as one that changes slowly over the
graph and of signals with large total variation as those that
change rapidly over the graph.
Consider the total variation of the eigenvectors vk and use
the facts that Lvk = λk vk and vkH vk = 1 to conclude that
TV(vk ) = vkH Lvk = λk .

(7)

It follows from (7) and the fact that the eigenvalues are
ordered as 0 = λ0 ≤ . . . ≤ λn−1 , that the total variations
of the eigenvectors vk follow the same order. Combining this
observation with the discussion following (6), we conclude
that when k is close to 0, the eigenvectors vk vary slowly
over the graph, whereas for k close to n the eigenvalues
vary more rapidly. Therefore, from (5) we see that the GFT
and iGFT allow us to decompose the brain signal x into
components that characterize different levels of variability.
The GFT coefficients x̃k for small values of k indicate how
much these slowly varying signals contribute to the observed
brain signal x. On the other hand, the GFT coefficients x̃k for
large values of k describe how much rapidly varying signals
contribute to the observed brain signal x.
B. Graph Filtering and Frequency Decompositions
Given a graph signal x with GFT x̃ we can isolate the
frequency components corresponding to the lowest KL graph
frequencies by defining the filtered spectrum x̃L = H̃L x̃
satisfying x̃Lk = x̃k for k < KL and x̃Lk = 0 otherwise. The
filter H̃L can be written as the diagonal matrix H̃L = diag(h̃L )
where the vector h̃L takes value 1 for frequencies smaller than
KL and is otherwise null,


h̃Lk = I k < KL .
(8)
Utilizing the definitions of the GFT in (3) and the iGFT in (4),
the spectral operation x̃L = H̃L x̃ is equivalent to performing
the following operations in the graph vertex domain
xL := Vx̃L = VH̃L x̃ = VH̃L V−1 x := HL x.

(9)

The last equality in (9) defines the matrix HL := VH̃L V−1 so
that the graph signal xL associated with low graph frequencies
of x can be written as the product xL = HL x. Since the signal
xL contains the low graph frequency components of x, we say
the matrix HL in (9) is a graph low-pass filter.
repreThe filter HL := VH̃L V−1 admits
Pn−1an alternative
k
sentation as the expansion HL =
h
L
in
terms
of
Lk
k=0
Laplacian powers [27]. The coefficients hLk in this expansion
are elements of the vector hL = Ψ−1 h̃L where Ψ is the
Vandermonde matrix defined by the eigenvalues of L, i.e.,


1
λ0
· · · λn−1
0

..
..  .
..
Ψ =  ...
(10)
.
.
. 
1

λn−1

···

λn−1
n−1

4

Since the eigenvalues are ordered in (10), the coefficients hLk
tend to P
be concentrated in small indexes k, and the expansion
n−1
HL = k=0 hLk Lk is therefore dominated by small powers
Lk . From this fact it follows that we can think of the graph
low-pass filtered signal xL as resulting from a localized averaging of the elements of x. To understand this interpretation,
simply note that L0 x = x coincides with the original signal,
Lx is an average of neighboring elements, L2 x is an average
of elements in nodes that interact via intermediate common
neighbors, and, in general, Lk x describes interactions between
k-hop neighbors. The fact that xL can be considered as a
signal that follows from local averaging of x implies that xL
has smaller total variation than x and is consistent with the
interpretation of low graph frequencies presented in Section
II-A. We point out that the definition hL = Ψ−1 h̃L assumes
the inverse matrix Ψ−1 exists. This holds true if the graph
Laplacian does not have duplicate eigenvalues, which is the
case for all functional brain networks examined in the paper.
Other types of graph filters can be defined analogously
to study interactions between signal components other than
the local interactions captured in xL . Apart from the graph
low-pass filter HL , we also consider a graph band-pass filter
HM and a graph high-pass filter HH , whose graph frequency
responses are defined as


h̃Mk = I KL ≤ k < KL + KM ,
(11)


h̃Hk = I KL + KM ≤ k .
(12)
The definitions in (8), (11), and (12) are such that the lowpass filter takes the lowest KL graph frequencies, the bandpass filter captures the middle KM graph frequencies, and the
high-pass filter the highest n − KL − KM frequencies. The
three filters are defined such that the graph frequencies of
their respective interest are mutually exclusive yet collectively
exhaustive. As a result, if we use xM := HM x and xH := HH x
to respectively denote the signals filtered by the band-pass and
high-pass filters, we have that the original signal can be written
as the sum x = xL + xM + xH . This gives a decomposition
of x into low, medium, and high frequency components
which respectively represent signals that have slow, medium,
and high variability with respect to the connectivity network
between brain regions. This decomposition is utilized in this
paper to analyze brain activity patterns associated with the
learning of visual-motor tasks.
III. B RAIN S IGNALS DURING L EARNING
We considered two experiments in which subjects learned
a simple motor task [28], [29]. In the experiments, fourtyseven right-handed participants (29 female, 18 male; mean age
24.13 years) volunteered with informed consent in accordance
with the University of California, Santa Barbara Internal
Review Board. After exclusions for task accuracy, incomplete
scans, and abnormal MRI, 38 participants were retained for
subsequent analysis.
Twenty individuals participated in the first experimental
framework. The experiment lasted 6 weeks, in which there
were 4 scanning sessions, roughly at the start of the experiment, at the end of the 2nd week, at the end of the

MIN Sequences
MOD Sequences
EXT Sequences

Session 1

Session 2

Session 3

Session 4

50
50
50

110
200
740

170
350
1430

230
500
2120

Fig. 1. Relationship between training duration, intensity, and depth for the
first experimental framework. The values in the table denote the number of
trials (i.e., “depth”) of each sequence type (i.e., “intensity”) completed after
each scanning session (i.e., “duration”) averaged over the 20 participants.

4th week, and at the end of the experiment, respectively.
During each scanning session, individuals performed a discrete
sequence-production task in which they responded to sequentially presented stimuli with their dominant hand on a custom
response box. Sequences were presented using a horizontal
array of 5 square stimuli with the responses mapped from
left to right such that the thumb corresponded to the leftmost
stimulus. The next square in the sequence was highlighted
immediately following each correct key press; the sequence
was paused awaiting the depression of the appropriate key
if an incorrect key was pressed. Each participant completed 6
different 10-element sequences. Each sequence consists of two
squares per key. Participants performed the same sequences at
home between each two adjacent scanning sessions, however,
with different levels of exposure for different sequence types.
Therefore, the number of trials completed by the participants
after the end of each scanning session depends on the sequence
type. There are 3 different sequence types (MIN, MOD, EXT)
with 2 sequences per type. The number of trials of each
sequence type completed after each scanning session averaged
over the 20 participants is summarized in Fig. 1. During
scanning sessions, each scan epoch involved 60 trials, 20 trials
for each sequence type. Each scanning session contained a
total of 300 trials (5 scan epochs) and a variable number of
brain scans depending on how quickly the task was performed
by the specific individual.
Eighteen subjects participated in the second experimental
framework. The experiment had 3 scanning sessions spanning
the three days. Each scanning session lasted roughly 2 hours
and no training was performed at home between adjacent
scanning sessions. Subjects responded to a visually cued
sequence by generating responses using the four fingers of
their nondominant hand on a custom response box. Visual
cues were presented as a series of musical notes on a pseudomusical staff with four lines such that the top line of the staff
mapped to the leftmost key pressed with the pinkie finger. Each
12-note sequence randomly ordered contained three notes per
line. Each training epoch involved 40 trials and lasted a total
of 245 repetition times (TRs), with a TR of 2,000 ms. Each
training session contained 6 scan epochs (240 trials) and lasted
a total of 2,070 scan TRs.
In both experiments participants were instructed to respond
promptly and accurately. Repetitions (e.g., “11”) and regularities such as trills (e.g., “121”) and runs (e.g., “123”) were
excluded in all sequences. The order and number of sequence
trials were identical for all participants. Participants completed
the tasks inside the MRI scanner for scanning sessions.
Reordering with fMRI was conducted using a 3.0 T Siemens

5

Fig. 3. Absolute magnitude at each of the n cortical structures averaged across participants in the 6 week experiment and averaged across all frequency
KL −1
KL +KM −1
components in (a) the set of low graph frequencies {vk }k=0
, (b) the set of middle graph frequencies {vk }k=K
, and (c) the set of high graph
L
n−1
frequencies {vk }k=K +K . (d)-(f) presents the average magnitudes for the 3 day experiment. Only brain regions with absolute magnitudes higher than a
L
M
fixed threshold (0.015) are colored. Magnitudes across the datasets are highly similar in the low and high graph frequencies (correlation coefficients 0.5818
and 0.6616, respectively). The brain regions with high magnitude values significantly overlap with the visual and sensorimotor modules, in which more than
60% of values greater than the threshold belong to the visual and sensorimotor modules.

Trio with a 12-channel phased-array head coil. For each
functional run, a single-shot echo planar imaging sequence
that is sensitive to blood oxygen level dependent (BOLD)
contrast was utilized to obtain 37 (the first experiment) or 33
(the second experiment) slices (3mm thickness) per repetition
time (TR), an echo time of 30 ms, a flip angle of 90◦ , a
field of view of 192 mm, and a 64 × 64 acquisition matrix.
Image preprocessing was performed using the Oxford Center
for Functional Magnetic Resonance Imaging of the Brain
(FMRIB) Software Library (FSL), and motion correction was
performed using FMRIB’s linear image registration tool. The
whole brain is parcellated into a set of n = 112 regions of
interest that correspond to the 112 cortical and subcortical
structures anatomically identified in FSL’s Harvard-Oxford
atlas. The choice of parcellation scheme is the topic of
several studies in structural [30], resting-state [31], and taskbased [32] network architecture. The question of the most
appropriate delineation of the brain into nodes of a network is
open and is guided by the particular question one wants to ask.
We use Harvard-Oxford atlas here because it is consistent with
previous studies of task-based functional connectivity during
learning [28], [29]. The threshold in probability cutoff settings
of Harvard Oxford atlas parcellation is 0 so that no voxels were
excluded.

coherence [33] between the activity of all possible pairs of
regions to construct n × n functional connectivity matrices
W. Besides, for each pair of brain regions i and j, we use tstatistical testing to evaluate the probability pi,j of observing
the measurements by random chance, when the actual data
are uncorrelated [34]. In the 3 day dataset, the value of all
elements with no statistical significance (pi,j > 0.05) [35] are
set to zero; the values remain unchanged otherwise. In the
3 day experiment, a single brain network is constructed for
each participant. Thresholding is applied because the networks
are for the entire span of the experiment and many entries
in W would be close to zero without threshold correction.
In the 6 week experiment, due to the long duration of the
experiment, we build a different brain network per scanning
session, per sequence type for each subject. Because each
network describes the functional connectivity for one training
session given a subject, not many entries will be removed
even in the presence of threshold correction; consequently, no
thresholding is applied for the 6 week dataset. We normalize
the regional mean BOLD observations x̂(t) at any sample time
t and consider x(t) = x̂(t)/kx̂(t)k2 such that the total energy
of activities at all structures is consistent at different t to avoid
extreme spikes due to head motion or drift artifacts in fMRI.

For each individual fMRI dataset, we estimate regional
mean BOLD time series by averaging voxel time series in each
of the n regions. We evaluate the magnitude squared spectral

IV. B RAIN N ETWORK F REQUENCIES
In this section, we analyze the graph spectrum brain networks of the dataset considered. For the brain network W

Weighted Zero Crossings

6

Total Variation

40
30
20
10
0

20

40

60

80

7000
6000
5000
4000
3000
2000
1000
0

100

20

40

(a)
Weighted Zero Crossings

Total Variation

15
10
5

40

60

100

300
250
200
150
100
50
0

20

80

(b)

20

0

60

Eigenvalue Index

Eigenvalue Index

80

100

20

40

60

80

100

Eigenvalue Index

Eigenvalue Index

(c)

(d)

Fig. 2. (a) Total variation TV(vk ) and (b) weighted zero crossings ZC(vk )
of the graph Laplacian eigenvectors for the brain networks averaged across
participants in the 6 week training experiment. (c) and (d) present the values
for the 3 day experiment. In both cases, the Laplacian eigenvectors associated
with larger indexes vary more on the network and cross zero relatively more
often, confirming the interpretation of the Laplacian eigenvalues as notions
of frequencies. Besides, note that total variation increases relatively linearly
with indexes.

of each subject, we construct its Laplacian L = D − W,
and evaluate the total variation TV(vk ) [cf. (7)] for each
eigenvector vk . Fig. 2 (a) and (c) plot the total variation of all
graph eigenvectors averaged across participants of the 6 week
training experiment and 3 day experiment, respectively. In both
experiments, the Laplacian eigenvectors associated with larger
indexes fluctuate more on the network. Another observation is
that the total variation increases in linear scale with indexes.
Besides total variation, the number of zero crossings is used
as a measure of the smoothness of signals with respect to an
underlying network [15]. Since brain networks are weighted,
we adapt a slightly modified version – weighted zero crossings
– to investigate the given graph eigenvector vk
1X
ZC(vk ) =
wij I {vk (i)vk (j) < 0} .
(13)
2
i6=j

In words, weighted zero crossings evaluate the weighted sum
of the set of edges connecting a vertex with a positive
signal to a vertex with a negative signal. Fig. 2 (b) and
(d) demonstrate the weighted zero crossings of all graph
eigenvectors averaged across subjects of the 6 week and 3
day experiments, respectively. The weighted zero crossings
increase proportionally with graph frequency index k for
0 ≤ k ≤ 100, as expected. However, eigenvectors associated
with higher graph frequencies exhibit lower weighted zero
crossings when k is greater than 100.
It would be interesting to examine where the associated
eigenvectors lie anatomically, and the relative strength of their
values. To facilitate the presentation, we consider three sets of
KL +KM −1
L −1
eigenvectors, {vk }K
and {vk }n−1
k=0 , {vk }k=KL
k=KL +KM ,
and compute the absolute magnitude at each of the n cortical
and subcortical regions averaged across participants and across

all graph frequencies belonging to each of the three sets. Fig.
3 presents the average magnitudes for the two experimental
frameworks considered in the paper using BrainNet [36],
where brain regions with absolute magnitudes lower than a
fixed threshold are not colored. Throughout the paper, the
parameter KL is set as 40 and KM is set as 32. We chose
this combination because it yields three roughly equally-sized
components with one piece corresponding to the 40 lowest
graph frequencies and another piece corresponding to the
40 highest frequencies. The results presented in the paper
are robust with the choice of parameters: we examined the
results for KL and KM in the range of 32 to 42 and found
similar observations as the ones presented. To demonstrate
that, we also consider the parameters KL = 38, KM = 37
and KL = 35, KM = 42. The average correlation coefficient
between the absolute magnitude across brain regions of the
parameters selected for this paper and the two chosen for
robustness is 0.8659. See Fig. 8 for another quantification
on the robustness of parameters. In this paper, we report
correlation coefficient as a quantification of similarity measure
when we examine the similarity between two vectors. We also
investigate cosine similarity and find high similarities as well.
A. Artificial Functional Brain Networks
An approach to analyze the complex networks is to define a model to generate artificial networks [37]. The main
motivation of an artificial network model is to use them to
analyze complex brain networks. Examples of such models
include Erdős-Rényi model [38] for unweighted networks,
Barabási-Albert model for scale free network [39], and recent
development and insights on weighted network models [37].
Here we present a framework to construct artificial networks
that can be used to mimic the functional brain networks with
only a few parameter inputs. The model is related to weighted
block stochastic model [40], but involves more aspects like
individual variance and links independent of their connecting
regions. The output of the method would be a symmetric
network with edge weights between 0 and 1 without self-loops.
To begin, suppose the desired network has two clusters of
nodes V1 and Vo . The algorithm requires the average edge
weight µ1 for connections between nodes of the first cluster
V1 , average edge weight µo for links between nodes of the
other cluster Vo , and average edge weight µ1o for intercluster connections. To reflect the fact that the edge weights
on some links are independent of their joining vertices, for
each edge within V1 , with probability p < 1, its weight
is randomly generated with respect to uniform distribution
U[0, 1] between 0 and 1, and with probability 1−p , its weight
is randomly generated with respect to uniform distribution
U[µ1 − δ, µ1 + δ]. The parameter p determines the percentage
of edges whose weights are selected irrespective of their actual
locations To further simulate the observation that different
participants may possess distinctive brain networks, if the
edge weight is randomly generated from a uniform distribution
U[µ1 −δ, µ1 +δ], it is then perturbed by wu ∼ U[−u /2, u /2]
where u controls the level of perturbation. The edge weights
for connections within cluster Uo are generated similarly:

7

with probability 1 − p , the edge weight is randomly chosen
from the uniform distribution U[µo − δ, µo + δ] before being
contaminated by wu ∼ U[−u /2, u /2]. The edge weights
for connections between clusters U1 and Uo are formed
analogously using µ1o . The method presented here can be
easily generalized to analyze brain networks with more regions
of interest, i.e. by specifying sets of regions of interest and
by detailing the expected correlation values on each type of
connection between different regions.
Remark 1 At one extreme we can make each node i belonging to a different set Vi = {i}. Then the method requires
the inputs of expected weights for all nodes, or alternatively
speaking, the expected network. At the other extreme, there
is only one set of nodes V, and then the method is highly
akin to a network with edge weights completely randomly
generated. Any construction of interest would have some prior
knowledge regarding the community structure. Therefore, the
method proposed here can be used to see if the network
constructed with the specific choice of community structure
highly simulate the key properties of the actual network, and
can be used to examine the evolution of community structure
in the brain throughout the process to master a particular task.
B. Spectral Properties of Brain Networks
In this section, we analyze graph spectral properties of brain
networks. Given the graph Laplacian, we examine the fluctuation of its eigenvectors on different types of connections in
the brain network [23]. More specifically, given an eigenvector
vk , its variation on the visual module is defined as
P
2
i,j∈Vv ,i6=j wij (vk (i) − vk (j))
P
,
(14)
TVvisual (vk ) =
i,j∈Vv ,i6=j wij
where Vv denotes the set of nodes belonging to the visual
module. The measure TVvisual (vk ) computes the difference for
signals on the visual module for each unit of edge weight.
To facilitate interpretation, we only consider three sets of
KL +KM −1
n−1
L −1
eigenvectors {vk }K
.
, and {vk }k=K
k=0 , {vk }k=KL
L +KM
We then compute the visual module total variation TVLvisual
M
L −1
averaged over eigenvectors {vk }K
k=0 , and TVvisual as well
H
as TVvisual similarly. Besides TVvisual , we also examine the
level of fluctuation of eigenvectors on edges within the motor
module, denoted as TVmotor , and on connections belonging to
brain modules other than the visual and motor module TVothers .
Further, there are links between two separate brain modules,
and to assess the variation of eigenvectors on those links, we
define total variations between the visual and motor modules
P
2
i∈Vv ,j∈Vm wij (vk (i) − vk (j))
P
TVvisual-motor (vk ) =
, (15)
i∈Vv ,j∈Vm wij
where Vm denotes the set of nodes belonging to the motor
module. Total variations TVvisual-others between the visual and
other modules, and total variations TVmotor-others between the
motor and other modules are defined analogously. We chose
to study visual and motor modules separately from other brain
modules because of their well-known associations with motor
learning [22], [23].

Fig. 4 (a) presents boxplots of the variation for eigenvectors
of different graph frequencies measured over different types
of connections across participants, at the start of the six week
training. Despite that total variation of eigenvectors should
increase with their frequencies, the variation on the other module TVLothers of eigenvectors associated with low frequencies
are higher than TVH
others (pass t-test with p < 0.0001). This
observation is discussed in detail in Section IV-C.
Next we study how the graph spectral properties of brain
networks evolve as participants become more familiar with
the tasks. Fig. 4 (b) illustrates the median of the variation
for eigenvectors of different graph frequencies measured over
different types of connections across subjects, at 10 different
levels of exposure in the six week training. As participants
become more acquainted with the assignment, their brain
networks display lower variation in the visual and motor
modules and higher variation in the other modules for low
and middle graph frequencies, and the exact opposite is true
for high graph frequencies. The association with training intensity is statistically significant (average correlation coefficient
r = 0.8164).
C. Discussion
Firstly, we examine why we see a decrease in zero crossings
of graph frequencies when k is greater than 100 in Fig. 2. A
detailed analysis shows this is because the functional brain
networks are highly connected with nearly homogeneous degree distribution, and consequently each high graph frequency
tends to have a value with high magnitude at one vertex of
high degree and similar values at other nodes, resulting in a
smaller global zero crossings for eigenvectors associated with
very high frequencies.
Secondly, in terms of the visualization of graph frequencies
in Fig. 3, the most interesting finding relates to the eigenvectors associated with high graph frequencies. The magnitudes
at different brain regions for high frequencies are significantly
similar across the two datasets investigated (correlation coefficient 0.6616), and brain regions with high magnitude values
are highly alike (greater than 60% overlap) to the visual and
sensorimotor cortices [41]. This is likely to be a consequence
of the fact that visual and motor regions are more strongly
connected with other structures, and hence an eigenvector with
a high magnitude on visual or motor structures would result
in high global spatial variation. The eigenvectors of low graph
frequencies are more spread across the networks, resulting
in low global variations. The middle graph frequencies are
less interesting – the magnitudes at most regions (greater
than 90%) do not pass the threshold, and little associations
(correlation coefficient 0.3529) can be found between the
eigenvectors of the 6 week and 3 day experiments.
Thirdly, to better interpret the meaning of variations for
specific types of connections, we construct artificial networks
as described in Section IV-A with visual and motor modules
as regions of interest, and consider other modules to be brain
regions other than visual and motor modules. We observe that
there are three contributing factors that cause the variation
within a specific module to become higher for higher eigenvectors and to become lower for lower eigenvectors: (i) Increases

Total Variation Magnitude

8

0.03

0.03

0.03

0.02

0.02

0.02

0.01

0.01

0.01

0

0
V

M

VM O
Region

Visual

Total Variation Magnitude

0

VO MO

V

Motor

M

Visual-Motor

Other

VO MO

V

Visual - Other

0.03

0.03

0.02

0.02

0.02

0.01

0.01

0.01

0

0

5

10

0

0
0

5

10

0

0.03

0.03

0.02

0.02

0.02

0.01

0.01

0.01

0

5
Trial

10

0

0
0

5
Trial

VM O
Region

VO MO

5

10

Trial

Trial

0.03

0

M

Motor - Other

0.03

Trial

Total Variation Magnitude

VM O
Region

10

0

5
Trial

10

Fig. 4. Spectral property of brain networks in the 6 week experiment. (a) Left: Averaged total variation of eigenvectors vk for 6 different types of connections
KL −1
of the brain averaged across all eigenvectors associated with low graph frequencies vk ∈ {vk }k=0
, across all participants and scan sessions. Middle: Across
KL +KM −1
all eigenvectors with middle graph frequencies vk ∈ {vk }k=K
. Right: Across all eigenvectors with high graph frequencies vk ∈ {vk }n−1
k=KL +KM . (b)
L
Median total variations of brain networks across participants of different scanning sessions and different sequence types with respect to the level of exposure
of participants to the sequence type at the scanning session. Relationship between training duration, intensity, and depth is summarized in Fig. 1. Value of
1 on the x-axis in the figures refers to minimum exposure to sequences (all 3 sequence types of the first session), and value of 10 on the x-axis denotes
the maximum exposure to sequences (EXT sequence types of the fourth session). An association between spectral property of brain networks and the level
of exposure is clearly observed (average correlation coefficient 0.8164). (c) Median total variations evaluated upon artificial networks. Spectral properties
of actual brain networks can be closely simulated using a few parameters. The main text gives all correlation values for similarity between variance among
subjects and between correlations of training intensity.

in the average edge weight for connections within the module,
(ii) Increments in the average edge weight for links between
this module to other module, and (iii) Escalation in the average
edge weight for associations within the other module. This
can also be observed by analyzing closely the definition of
total variation. If a module is highly connected, in order for
the eigenvector associated with a low graph frequency to be
smooth on the entire network, it has to be smooth on the
specific module, resulting in a low value in the variation of
an eigenvector associated with a low graph frequency with
respect to the module of interest. Similarly, the increase in the

variation of connections between two modules, e.g. between
visual and other modules are resulted from: (i) The growth in
the average edge weight for connections between visual and
other modules, or (ii) The augmentation of average weight for
links within the other module. The graph spectral properties as
in Fig. 4 (a) are observed because (i) visual and motor modules
are themselves highly connected, and (ii) visual module is also
strongly linked with motor module.
Finally, in analyzing the evolution of graph spectral properties as participants become more familiar with the tasks, following the interpretations based on artificial network analysis,

9

Fig. 5. Distribution of decomposed signals for the 6 week experiment. (a) Absolute magnitudes for all brain regions with respect to xL – brain signals
varing smoothly across the network – averaged across all sample points for each individual and across all participants at the first scan session of the 6 week
dataset. (b) With respect to xM and (c) with respect to xH – signals rapidly fluctuating across the brain. (d), (e), and (f) are averaged xL , xM and xH at the
last scan session of the 6 week dataset, respectively. Only regions with absolute magnitudes higher than a fixed threshold are colored.

Fig. 6. Distribution of decomposed signals for the 3 day experiment. (a), (b), and (c) are xL , xM and xH averaged across all sample points for each subject
and across participants in the 3 day experiment, respectively. Regions with absolute value less than a threshold are not colored.

this evolution in graph spectral properties of brain networks
is mainly caused by the decrease in values of connections
within visual and motor modules and between the visual and
motor modules. An interesting observation is that the values in
the variation of eigenvectors associated with high frequencies
decline with respect to the visual module much faster than
that of motor module, even though the visual module is more
strongly connected throughout training compared to the motor
module. A deep analysis using artificial networks shows that
this results from the following three factors: (i) Though more
strongly connected compared to the motor module, connections within the visual module weaken very quickly, (ii) The

motor module is more closely connected with the other module
than the link between the visual module to the other module,
and (iii) Association levels within the other module stay
relatively constant. Therefore, as participants become more
exposed to the tasks, compared to the visual module, the motor
module becomes more strongly connected. The graph spectral
properties of actual brain networks and their evolution can be
closely imitated using artificial networks as plotted in Fig. 4
(c). The artificial network created for our analysis best imitated
the real brain networks with parameters p of 0.10, u of 0.10,
and δ of 0.01. The average edge weights µ for visual (v), motor
(m), other (o), and inter-connecting regions are µv = 0.6028,

10

Fig. 7. Temporal adaptations of spatial variations. Boxplots showing differences in temporal adaptabilities between brain activities with smooth (pink),
moderate (red) and rapid (maroon) spatial variations, measured over the complete experiment (Left) and individual training sessions (Right) for 6 week
experiment (Top) and 3 day experiment (Bottom). We measured the temporal adaptations using the variance of the averaged activities over the complete
experiment or with individual training sessions. Compared to activities with moderate spatial variations, smooth (95% sessions pass t-test with p < 0.01) and
rapid (65% sessions pass t-test with p < 0.005) spatial variations have significantly higher temporal adaptations.

µm = 0.4902, µo = 0.3098, µvm = 0.3985, µvo = 0.3181,
and µmo = 0.3271. The correlation coefficients of association
with training intensity between real and artificial networks
for low, medium, and high graph frequencies are 0.6436,
0.7187 and 0.8457, respectively. Additionally, the variation
among participants in real dataset can be highly mimicked
using artificial network model we proposed, with correlation
coefficients 0.9338, 0.9660, and 0.9486 for low, medium, and
high graph frequencies, respectively. The analysis for the three
day training dataset is highly similar (correlation coefficients
0.9834, 0.9186, and 0.9674 for low, medium, and high graph
frequencies, respectively) and for this reason we do not present
and analyze it separately here.

V. F REQUENCY D ECOMPOSITION OF B RAIN S IGNALS
The previous sections focus on the study of brain networks
and their graph spectral properties. In this section, we investigate brain signals from a GSP perspective, and analyze
the brain signals by examining the decomposed graph signals
xL , xM , and xH with respect to the underlying brain networks.
We compute the absolute magnitude of the decomposed signal
xL for each brain region averaged across all sample signals for
each individual during a scan session and then averaged across
all participants. Similar aggregation is applied for xM and xH .
Fig. 5 presents the distribution of the decomposed signals
corresponding to different levels of spatial variations for the
first scan session (top row) and the last scan session (bottom
row) in the 6 week experiment. Fig. 6 exhibits how the
decomposed signals are distributed across brain regions in
the 3 day experiment. Brain regions with absolute magnitudes
lower than a fixed threshold are not colored.

A. Temporal Variation of Graph Frequency Components
We analyze temporal variation of decomposed signals with
respect to different levels of spatial variations. To that end, we
evaluate the variance of the decomposed signals over multiple
temporal scales – over days and minutes [42] – for the two
experiments. We describe the method specifically for xL for
simplicity and similar computations were conducted for xM
and xH . At the macro timescale, we average the decomposed
signals xL for all sample points within each scanning session
with different sequence type, and evaluate the variance of
the magnitudes of the signals [16] across all the scanning
sessions and sequence types. For the 6 week experiment,
there are 4 scanning sessions and 3 different sequence types,
so the variance is with respect to 12 points. For the 3 day
experiment, there are 3 scanning sessions and only 1 sequence
type, and therefore the variance is for 3 points. As for the
micro or minute-scale, we average the decomposed signals
xL for all sample points within each minute, and evaluate
the variance of the magnitudes of the averaged signals across
all minute windows for each scanning session with different
sequence types. The evaluated variance is then averaged across
all participants of the experiment of interest.
Fig. 7 displays the variance of the decomposed signals
xL , xM and xH at two different temporal scales of the two
experiments. For the 6 week dataset, 3 session-sequence
combinations, with the number proportional to the level of
exposure of participants to the sequence (1-MIN refers to MIN
sequence at session 1, 5 denotes MIN sequence at session 4, 9
entails EXT sequence at session 3) are selected out of the 12
combinations in total for a cleaner illustration, but all the other
session-sequence combinations exhibit similar properties.

11

B. Discussion
A deep analysis of Figs. 5 and 6 yields many interesting
aspects of graph frequency decomposition. First, for xL , the
magnitudes on adjacent brain regions tend to possess highly
similar values, resulting in a more evenly spread brain signal
distribution, where as for xH , neighboring signals can exhibit
highly dissimilar values; this corroborates the motivation to
use graph frequency decomposition to segment brain signals
into pieces corresponding to different levels of spatial fluctuations. Second, decomposed signal of a specific level of
variations, e.g. xL and xH , are highly similar with respect
to different scan sessions in an experiment as well as with
respect to two experiments with different sets of participants.
The correlation coefficient between datasets for low, medium
and high graph frequencies are 0.5841, 0.2852, and 0.6469, respectively. This reflects the fact that frequency decomposition
is formed by applying graph filters with different pass bands
upon signals and therefore should express some consistent
aspects of brain signals. Third, recall that we normalize the
brain signals at every sample point for all subjects, and for this
reason signals xL , xM and xH would be similarly distributed
across the brain if nothing interesting happens at the decomposition. However, in both Fig.s 5 and 6, it is observed that
many brain regions possess magnitudes higher than a threshold
in xL (∼ 60% pass) and xH (∼ 20% pass) while not many
brain regions pass the thresholding with respect to xM (∼ 3%
pass). It has long been understood that the brain combines
some degree of disorganized behavior with some degree of
regularity and that the complexity of a system is high when
order and disorder coexist [43]. xL varies smoothly across
the brain network and therefore can be regarded as regularity
(order), whereas xH fluctuates rapidly and consequently can be
considered as randomness (disorder). This evokes the intuition
that graph frequency decomposition segments a brain signal x
into pieces xL and xH , which reflect order and disorder (and
are therefore more interesting), as well as the remaining xM .
For the variance analysis, it is expected for the low graph
frequency components (smooth spatial variation) to exhibit
the smallest temporal variations, exceeded by medium and
then high counterparts. Nonetheless, it is observed that brain
activities with smooth spatial variations exhibit the most
rapid temporal variation. Because it has been shown that
temporal variation of observed brain activities is associated
with better performance in tasks [16], this indicates a stronger
contribution of low graph frequency components during the
learning process. Furthermore, since the measurements were
normalized such that the total energy of overall brain activities
stayed constant at different sampling points, the rapid temporal
changes of low graph frequency components should be accompanied by fast temporal variation of some other components,
which are found to be high frequency components in all cases.
Because these results were consistent for all of the temporal
scales and datasets that we examined, and the association
between temporal variability and positive performance has
been established [17], we concluded that brain activities with
smooth or rapid spatial variations offer greater contributions
during learning. The graph frequency signatures at different

6 week experiment (linear scale)
6 week experiment (logarithm scale)
3 day experiment

kxL k2

kxM k2

kxH k2

−0.3155
−0.5409
−0.9873

0.0897
0.3992
0.8443

0.4125
0.3565
0.9605

Fig. 8. Pearson correlation coefficients between the number of trials (level
of task familiarity) and R values, defined as correlations between learning
rate parameters and the norm of the decomposed signal of interest. More
obvious adaptability between decomposed signals and learning across training
is observed for xL and xH , with decreasing association with exposure to tasks
for the former and increasing importance for the latter.

stages of learning is analyzed in the next section.
VI. F REQUENCY S IGNATURES OF TASK FAMILIARITY
Given that the decomposed signals exhibit interesting perspectives, it is natural to probe whether the signals corresponding to different levels of spatial variations associate with
learning. To that end, we first describe how learning rate is
evaluated. Given a participant, for each sequence completed,
we defined the movement time M as the difference between
the time of the first button press and the time of the last
button press during a single sequence. We then estimate the
participant’s learning rate by fitting an exponential function
(plus a constant) using the robust outlier correction [44] to
the sequence of movement times M
M = c1 et/κ + c2 .

(16)

where t is a sequence representing the time index, κ is the
exponential drop-off parameter (which we call the “learning
rate parameter”) used to describe the early and fast rate of
improvement, and c1 and c2 are nonnegative constants. Their
sum c1 + c2 is an estimation of the starting speed of the
participant of interest prior to training, while the parameter
c2 entails the fastest speed to complete the sequence attained
by that participant after extended training. A negative value
of κ indicates a decrease in movement time M (t), which
is thought to indicate that learning is occurring [45]. We
chose exponential because it is viewed as the most statistically
robust choice [46]. Further, the approach that we used has the
advantage of estimating the rate of learning independent of
initial performance or performance ceiling.
We evaluate the learning rate for all participants at each
scanning session, and then compute the correlation between
the norm kxL k2 of the decomposed signal corresponding to
low spatial variation and the learning rates across subjects.
The correlation (R value) between the norms kxM k2 as well
as kxH k2 and learning rates are also calculated. Fig. 9 plots
the Pearson correlation coefficients at all scanning sessions of
the two experiments considered. The horizontal axis denotes
the level of exposure of participants to the sequence – which
day in the 3 day experiment and how many number of trials
participants have completed at the end of the scanning session
in the 6 week experiment. Points are densely distributed for
small number of trials in the 6 week experiment, so to mitigate
this effect, we also plot the points by taking the logarithm
of numbers of trials completed. We emphasize that due to
normalization at each sampling point, the correlation values

12

0.5

0.5

0

0

0

-0.5

-0.5

-0.5

R Value

0.5

0

500 1000 1500 2000 2500

0

500 1000 1500 2000 2500

R Value

Number of Trials
0.5

0.5

0

0

0

-0.5

25

100

1000

-0.5

25

100

Number of Trials

1000

25

0.5

0

0

0

-0.5

-0.5
3

1000

Number of Trials

0.5

2

100

Number of Trials

0.5

1

500 1000 1500 2000 2500

Number of Trials

0.5

-0.5

R Value

0

Number of Trials

-0.5
1

Day

2

3

Day

1

2

3

Day

Fig. 9. Scatter plots depicting the number of trials (level of task familiarity) and R values, defined here as correlations between learning rate parameters
and the norm of the decomposed signal of interest (Pink points in the Left: xL , Red points in the Middle: xM , and Maroon points in the Right: xH ). Top
row: 6 week experiment with number of trials described in linear scale. Middle row: 6 week experiment withe number of trials evaluated in logarithm scale.
We examine 6 week experiment by ordering the number of trials in both linear and logarithm scales to alleviate the fact that number of trials are densely
distributed towards small values. Bottom row: 3 day experiment with scanning session ordered by the number of days in the experiment.

would all be 0 if graph frequency decomposition segments
brain signals into three equivalent pieces. There are scan
sessions where the correlation is of particular interest, however
the most noteworthy observation is the change of correlation
values with the level of exposure of participants. For xL
corresponding to smooth spatial variation, its correlation with
learning is above zero (≈ 0.25) at the start of the training
when participants perform the task for the first time. The
correlation gradually decreases with heavier training intensity
until below zero (≈ −0.25) at the end of the experiment
when individuals are highly familiar with the sequence. For xH
corresponding to vibrant spatial variation, its correlation with
learning is below zero (≈ −0.2) at the start of the training, and
gradually increases throughout training until it is above zero
(≈ 0.25) at the end of the experiment – the exact opposite of
xL . For xM , correlation between its norm kxM k2 with learning
rate generally increases with the intensity of training, however
this trend is not as obvious compared to other decomposition
counterparts. The correlation between the number of trials and
R value is summarized in Fig. 8. For robustness testing, we
conduct similar analysis using the other two sets of parameters
described in Section IV. The average absolute difference for
entries in Fig. 8 between the original parameters and the other
two parameter choices is 0.0377. Again, similar observations
are found in different experiments involving different learning
tasks and different sets of participants.

A. Discussion
This result further implies that the most association between
learning or adaptability during the training process comes from
the brain signals that either vary smoothly (xL , regularity) or
rapidly (xH , randomness) with respect to the brain network.
Therefore, the graph frequency decomposition could be used
to capture more informative brain signals by filtering out noninformative counterparts, most likely associated with middle
graph frequencies. Besides, the positive association between
kxL k2 and learning rates as well as the negative association
between kxH k2 and learning rates at the start of training
indicates that it favors learning to have more smooth, spread,
and cooperative brain signals when we face an unfamiliar
task. As we gradually become familiar with the task, the
smooth and cooperative signal distribution becomes less and
less important, and there is a level of exposure when such
signal distribution becomes destructive instead of constructive.
We note that the task in the 3 day experiment is more difficult
compared to that of the 6 week experiment, and therefore the
time when the cooperative signal distribution starts to become
detrimental (the point where the regression line intercepts the
horizontal line of R value equaling 0) is also comparable in
the two experiments, describing a certain level of familiarity
to the task. When we become highly familiar with the task,
it is better and favors further learning to have varied, spiking,
and competitive brain signals.
In the dataset evaluated here, we utilize the average co-

13

herence between time series at pairs of brain cortical and
subcortical regions during the training as the network. Hence,
a concentration of brain activities towards low graph frequencies would imply that activities on brain regions that
are generally cooperative are indeed similar. Simultaneously,
the interpretation of concentration of brain activities towards
high graph frequencies is that brain activities on brain regions
that are generally cooperative are in fact dissimilar. In terms
of learning, one possible explanation is that there are two
different stages in learning: we start by grasping the big picture
of the task to perform relatively well, and then we refine the
details to perform better and to approach our limits.
Because the graph frequency analysis method presented in
this paper applies to any setting where signals are defined on
top of a network structure representing proximities between
nodes, it would be interesting in future to use this method to
investigate other types of signals and networks in neuroscience
problems. As an example, in situations given fMRI measurements on structural networks, concentration of signals in low
graph frequency components would imply functional activities
do behave according to the structural networks.
Besides, it has been understood that learning is different
when one is unfamiliar or familiar with a particular task – it
is easy to improve performance at first exposure due to the fact
that one is far from their performance ceiling. It would therefore be interesting to utilize graph frequency decomposition
to further analyze the difference between learning scenarios at
different stages of familiarity, e.g. adaptability at first exposure
and creativity when one fully understands the components of
the specific tasks.
VII. C ONCLUSION
We used graph spectrum methods to analyze functional
brain networks and signals during simple motor learning tasks,
and established connections between graph frequency with
principal component analysis when the networks of interest
denote functional connectivity. We discerned that brain activities corresponding to different graph frequencies exhibit
different levels of adaptability during learning. Further, the
strong correlation between graph spectral property of brain
networks with the level of familiarity of tasks was observed,
and the most contributing frequency signatures at different task
familiarity was recognized.
R EFERENCES
[1] H. Haken, Principles of brain functioning: a synergetic approach to
brain activity, behavior and cognition. Springer Science & Business
Media, 2013, vol. 67.
[2] M. D. Fox and M. E. Raichle, “Spontaneous fluctuations in brain activity
observed with functional magnetic resonance imaging,” Nature Reviews
Neuroscience, vol. 8, no. 9, pp. 700–711, 2007.
[3] J. Britz, D. Van De Ville, and C. M. Michel, “Bold correlates of eeg
topography reveal rapid resting-state network dynamics,” Neuroimage,
vol. 52, no. 4, pp. 1162–1170, 2010.
[4] L. J. Cole, M. J. Farrell, E. P. Duff, J. B. Barber, G. F. Egan, and
S. J. Gibson, “Pain sensitivity and fmri pain-related brain activity in
alzheimer’s disease,” Brain, vol. 129, no. 11, pp. 2957–2965, 2006.
[5] J. D. Medaglia, W. Huang, S. Segarra, C. Olm, J. Gee, M. Grossman,
A. Ribeiro, C. T. McMillan, and D. S. Bassett, “Brain network efficiency is influenced by pathological source of corticobasal syndrome,”
Neurology, vol. (submitted), 2016.

[6] S. Achard, R. Salvador, B. Whitcher, J. Suckling, and E. Bullmore, “A
resilient, low-frequency, small-world human brain functional network
with highly connected association cortical hubs,” The Journal of neuroscience, vol. 26, no. 1, pp. 63–72, 2006.
[7] E. Bullmore and O. Sporns, “The economy of brain network organization,” Nature Reviews Neuroscience, vol. 13, no. 5, pp. 336–349, 2012.
[8] J. Richiardi, S. Achard, H. Bunke, and D. Van De Ville, “Machine learning with brain graphs: predictive modeling approaches for functional
imaging in systems neuroscience,” Signal Processing Magazine, IEEE,
vol. 30, no. 3, pp. 58–70, 2013.
[9] H. Q. Nguyen, P. Chou, Y. Chen et al., “Compression of human body
sequences using graph wavelet filter banks,” in Acoustics, Speech and
Signal Processing (ICASSP), 2014 IEEE International Conference on.
IEEE, 2014, pp. 6152–6156.
[10] J. Ma, W. Huang, S. Segarra, and A. Ribeiro, “Diffusion filtering for
graph signals and its use in recommendation systems,” in Proc. Int. Conf.
Acoustics Speech Signal Process, vol. (submitted), Shanghai, China,
March 20-25 2016.
[11] S. Segarra, W. Huang, and A. Ribeiro, “Diffusion and superposition
distances for signals supported on networks,” IEEE Trans. Signal Inform.
Process. Networks, vol. 1, no. 1, pp. 20–32, March 2015.
[12] A. Gadde, A. Anis, and A. Ortega, “Active semi-supervised learning
using sampling theory for graph signals,” in Proceedings of the 20th
ACM SIGKDD international conference on Knowledge discovery and
data mining. ACM, 2014, pp. 492–501.
[13] A. Sandryhaila and J. M. Moura, “Discrete signal processing on graphs,”
IEEE Trans Signal Process, vol. 61, no. 7, pp. 1644–1656, 2013.
[14] ——, “Discrete signal processing on graphs: Frequency analysis,” IEEE
Trans Signal Process, vol. 62, no. 12, pp. 3042–3054, 2014.
[15] D. Shuman, S. K. Narang, P. Frossard, A. Ortega, P. Vandergheynst
et al., “The emerging field of signal processing on graphs: Extending
high-dimensional data analysis to networks and other irregular domains,”
IEEE Signal Process Mag, vol. 30, no. 3, pp. 83–98, 2013.
[16] D. D. Garrett, N. Kovacevic, A. R. McIntosh, and C. L. Grady, “The
modulation of bold variability between cognitive states varies by age
and processing speed,” Cerebral Cortex, p. bhs055, 2012.
[17] J. J. Heisz, J. M. Shedden, and A. R. McIntosh, “Relating brain signal
variability to knowledge representation,” Neuroimage, vol. 63, no. 3, pp.
1384–1392, 2012.
[18] N. Leonardi, J. Richiardi, M. Gschwind, S. Simioni, J.-M. Annoni,
M. Schluep, P. Vuilleumier, and D. Van De Ville, “Principal components
of functional connectivity: a new approach to study dynamic brain
connectivity during rest,” NeuroImage, vol. 83, pp. 937–950, 2013.
[19] R. Viviani, G. Grön, and M. Spitzer, “Functional principal component
analysis of fmri data,” Human brain mapping, vol. 24, no. 2, pp. 109–
129, 2005.
[20] L. M. Harrison, W. Penny, J. Daunizeau, and K. J. Friston, “Diffusionbased spatial priors for functional magnetic resonance images,” Neuroimage, vol. 41, no. 2, pp. 408–423, 2008.
[21] J. Shi and J. Malik, “Normalized cuts and image segmentation,” Pattern
Analysis and Machine Intelligence, IEEE Transactions on, vol. 22, no. 8,
pp. 888–905, 2000.
[22] J. A. Kleim, S. Barbay, N. R. Cooper, T. M. Hogg, C. N. Reidel, M. S.
Remple, and R. J. Nudo, “Motor learning-dependent synaptogenesis is
localized to functionally reorganized motor cortex,” Neurobiology of
learning and memory, vol. 77, no. 1, pp. 63–77, 2002.
[23] D. S. Bassett, M. Yang, N. F. Wymbs, and S. T. Grafton, “Learninginduced autonomy of sensorimotor systems,” Nature neuroscience,
vol. 18, no. 5, pp. 744–751, 2015.
[24] F. Chung, Spectral graph theory. American Mathematical Soc., 1997,
vol. 92.
[25] W.-T. Zhang, Z. Jin, G.-H. Cui, K.-L. Zhang, L. Zhang, Y.-W. Zeng,
F. Luo, A. C. Chen, and J.-S. Han, “Relations between brain network
activation and analgesic effect induced by low vs. high frequency electrical acupoint stimulation in different subjects: a functional magnetic
resonance imaging study,” Brain research, vol. 982, no. 2, pp. 168–178,
2003.
[26] J. Grutzendler, N. Kasthuri, and W.-B. Gan, “Long-term dendritic spine
stability in the adult cortex,” Nature, vol. 420, no. 6917, pp. 812–816,
2002.
[27] S. Segarra, A. G. Marques, G. Leus, and A. Ribeiro, “Reconstruction
of graph signals through percolation from seeding nodes,” Trans.
Signal. Process., vol. (submitted), 2015. [Online]. Available: http:
//arxiv.org/abs/1507.08364
[28] D. S. Bassett, N. F. Wymbs, M. P. Rombach, M. A. Porter, P. J. Mucha,
and S. T. Grafton, “Task-based core-periphery organization of human
brain dynamics,” PLoS Comput. Biol, vol. 9, no. 9, p. e1003171, 2013.

14

[29] D. S. Bassett, N. F. Wymbs, M. A. Porter, P. J. Mucha, J. M. Carlson,
and S. T. Grafton, “Dynamic reconfiguration of human brain networks
during learning,” Proc Amer Philos Soc Natl Acad Sci USA, vol. 108,
no. 18, pp. 7641–7646, 2011.
[30] A. Zalesky, A. Fornito, I. H. Harding, L. Cocchi, M. Yücel, C. Pantelis,
and E. T. Bullmore, “Whole-brain anatomical networks: does the choice
of nodes matter?” Neuroimage, vol. 50, no. 3, pp. 970–983, 2010.
[31] J. Wang, L. Wang, Y. Zang, H. Yang, H. Tang, Q. Gong, Z. Chen,
C. Zhu, and Y. He, “Parcellation-dependent small-world brain functional
networks: A resting-state fmri study,” Hum Brain Mapp, vol. 30, no. 5,
pp. 1511–1523, 2009.
[32] J. D. Power, A. L. Cohen, S. M. Nelson, G. S. Wig, K. A. Barnes, J. A.
Church, A. C. Vogel, T. O. Laumann, F. M. Miezin, B. L. Schlaggar
et al., “Functional network organization of the human brain,” Neuron,
vol. 72, no. 4, pp. 665–678, 2011.
[33] F. T. Sun, L. M. Miller, and M. D’Esposito, “Measuring interregional
functional connectivity using coherence and partial coherence analyses
of fmri data,” Neuroimage, vol. 21, no. 2, pp. 647–658, 2004.
[34] Y. He, Z. J. Chen, and A. C. Evans, “Small-world anatomical networks in
the human brain revealed by cortical thickness from mri,” Cereb cortex,
vol. 17, no. 10, pp. 2407–2419, 2007.
[35] C. R. Genovese, N. A. Lazar, and T. Nichols, “Thresholding of statistical maps in functional neuroimaging using the false discovery rate,”
Neuroimage, vol. 15, no. 4, pp. 870 – 878, 2002.
[36] M. Xia, J. Wang, Y. He et al., “Brainnet viewer: a network visualization
tool for human brain connectomics,” PloS one, vol. 8, no. 7, p. e68910,
2013.
[37] C. Lohse, D. S. Bassett, K. O. Lim, and J. M. Carlson, “Resolving
anatomical and functional structure in human brain organization: Identifying mesoscale organization in weighted network representations,”
2014.
[38] P. ERDdS and A. R&WI, “On random graphs i,” Publ. Math. Debrecen,
vol. 6, pp. 290–297, 1959.
[39] R. Albert and A.-L. Barabási, “Statistical mechanics of complex networks,” Reviews of modern physics, vol. 74, no. 1, p. 47, 2002.
[40] C. Aicher, A. Z. Jacobs, and A. Clauset, “Learning latent block structure
in weighted networks,” Journal of Complex Networks, p. cnu026, 2014.
[41] V. Braitenberg and A. Schüz, Cortex: statistics and geometry of neuronal
connectivity. Springer Science & Business Media, 2013.
[42] K. M. Newell, G. Mayer-Kress, S. L. Hong, and Y.-T. Liu, “Adaptation
and learning: Characteristic time scales of performance dynamics,” Hum
Mov Sci, vol. 28, no. 6, pp. 655–687, 2009.
[43] O. Sporns, Networks of the Brain. MIT press, 2011.
[44] D. A. Rosenbaum, Human motor control. Academic press, 2009.
[45] E. Dayan and L. G. Cohen, “Neuroplasticity subserving motor skill
learning,” Neuron, vol. 72, no. 3, pp. 443–454, 2011.
[46] A. Heathcote, S. Brown, and D. Mewhort, “The power law repealed:
The case for an exponential law of practice,” Psychonomic bulletin &
review, vol. 7, no. 2, pp. 185–207, 2000.

