NeuroImage 133 (2016) 251–265

Contents lists available at ScienceDirect

NeuroImage
journal homepage: www.elsevier.com/locate/ynimg

Estimating repetitive spatiotemporal patterns from resting-state brain
activity data
Yusuke Takeda ⁎, Nobuo Hiroe, Okito Yamashita, Masa-aki Sato
Department of Computational Brain Imaging, ATR Neural Information Analysis Laboratories, 2-2-2 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-0288, Japana r t i c l e

i n f oArticle history:
Received 31 July 2015
Accepted 4 March 2016
Available online 12 March 2016
Keywords:
Spatiotemporal pattern
Resting-state
Spontaneous
MEG
EEG
fMRI

a b s t r a c t
Repetitive spatiotemporal patterns in spontaneous brain activities have been widely examined in non-humanstudies. These studies have reported that such patterns reﬂect past experiences embedded in neural circuits. In
human magnetoencephalography (MEG) and electroencephalography (EEG) studies, however, spatiotemporal
patterns in resting-state brain activities have not been extensively examined. This is because estimating spatiotemporal patterns from resting-state MEG/EEG data is difﬁcult due to their unknown onsets. Here, we propose amethod to estimate repetitive spatiotemporal patterns from resting-state brain activity data, including MEG/EEG.
Without the information of onsets, the proposed method can estimate several spatiotemporal patterns, even if
they are overlapping. We veriﬁed the performance of the method by detailed simulation tests. Furthermore,
we examined whether the proposed method could estimate the visual evoked magnetic ﬁelds (VEFs) without
using stimulus onset information. The proposed method successfully detected the stimulus onsets and estimatedthe VEFs, implying the applicability of this method to real MEG data. The proposed method was applied toresting-state functional magnetic resonance imaging (fMRI) data and MEG data. The results revealed informative
spatiotemporal patterns representing consecutive brain activities that dynamically change with time. Using thismethod, it is possible to reveal discrete events spontaneously occurring in our brains, such as memory retrieval.
© 2016 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license(http://creativecommons.org/licenses/by/4.0/).

Introduction
Over the past decade, resting-state (or spontaneous) brain activities
have attracted much interest, prompting a growing body of neurosciencestudies.
In non-human studies, repetitive spatiotemporal patterns emerging
in spontaneous brain activities have been widely examined (e.g. Ikegaya
et al., 2004). Here, we deﬁne spatiotemporal patterns as activities
represented by two-dimensional matrices of channel × time (Fig. 1A).
These studies have reported that the repetitive spatiotemporal patterns
resemble the preceding brain activities during tasks, suggesting thatthese patterns reﬂect past experiences embedded in neural circuits(Foster and Wilson, 2006; Han et al., 2008; Ji and Wilson, 2007;
Wilson and McNaughton, 1994). Furthermore, it has been reported
that the spatiotemporal patterns are predictive of future brain activities
during tasks, suggesting that they contribute to the encoding of future
novel experiences (Dragoi and Tonegawa, 2011, 2013). Theoreticalstudies implied that encoding information by spatiotemporal patterns
has advantages in terms of the computational efﬁciency of patternrecognition and memory capacity (Hopﬁeld, 1995; Izhikevich, 2006).

⁎ Corresponding author. Fax: +81 774 95 1259.
E-mail address: takeda@atr.jp (Y. Takeda).

All of the above studies highlight the signiﬁcance of examining spatiotemporal patterns in spontaneous brain activities.
In human studies, functional connectivities of resting-state functional magnetic resonance imaging (fMRI) have been widely examined(Beckmann et al., 2005; Biswal et al., 1995, 2010; Fox et al., 2005; Fox
and Raichle, 2007; Raichle et al., 2001; Smith et al., 2009). Functional
connectivity is the correlation of fMRI time series across regions, and it
is examined using either seed-based correlation (Biswal et al., 2010)
or independent component analysis (ICA) (Beckmann et al., 2005).
Starting with motor cortices (Biswal et al., 1995), several sets of correlated brain regions were identiﬁed, such as the default mode network(DMN) (Fox et al., 2005; Raichle et al., 2001). Recently, Majeed et al.
(2011) developed a template-matching algorithm to estimate
spatiotemporal patterns from resting-state fMRI data, and they revealed
spatiotemporal patterns consisting of an alteration between activationof areas comprising the DMN and the task-positive network.
In human magnetoencephalography (MEG) and electroencephalography (EEG) studies, however, spatiotemporal patterns in resting-state
brain activities have not been thoroughly examined. Many of theresting-state MEG/EEG studies also focused on functional connectivity,
which here is the correlation of an oscillation's amplitudes acrossregions. The functional connectivities of MEG/EEG were reported to relate to those of fMRI (Baker et al., 2014; Brookes et al., 2011; de Pasquale
et al., 2010, 2012; Mantini et al., 2007). Using clustering methods, some

http://dx.doi.org/10.1016/j.neuroimage.2016.03.014
1053-8119/© 2016 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

252

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

studies examined spatial patterns in resting-state EEG data, the socalled microstates (Britz et al., 2010; Van de Ville et al., 2010). However,they did not examine the spatiotemporal patterns.
This paucity of research efforts devoted to examining the spatiotemporal patterns in resting-state MEG/EEG data could be attributed to the
difﬁculty of estimating them. In the case of task-related data, spatiotemporal patterns are conventionally estimated using externally observable
onsets, such as stimulus and response onsets. For example, thespatiotemporal MEG patterns time-locked at visual stimuli, that is, theso-called visual evoked magnetic ﬁelds (VEFs), are obtained by averaging MEG data triggered at the visual stimulus onsets. In the case ofresting-state data, however, externally observable onsets do not exist,
so the averaging procedure cannot be used. Considering that taskrelated MEG/EEG data exhibit spatiotemporal patterns time-locked toovert events, it is reasonable to assume that resting-state MEG/EEG
data also exhibit spatiotemporal patterns time-locked to covert events,
such as memory retrieval (Deuker et al., 2013; Staresina et al., 2013;
Tambini et al., 2010).
In the above method of Majeed et al. (2011), a segment starting at a
random time point is regarded as a template. Then, segments
resembling the template are searched for over time, and the templateis updated by averaging the found segments. This method works wellwhen the spatiotemporal patterns do not overlap. If they do overlap,
however, the method estimates the spatiotemporal patterns that arecontaminated by each other. In the cases of MEG/EEG data, signals
from different brain regions are spatially mixed, and thus spatiotemporal patterns probably do overlap. Therefore, their method is not suitablefor MEG/EEG data. Because MEG/EEG can capture brain activities on areal electrical time scale, estimating spatiotemporal patterns fromthese data is also important.
In this study, we propose a method to estimate repetitive spatiotemporal patterns from resting-state brain activity data, including MEG/

A

EEG. Without the information of onsets, the proposed method canestimate several spatiotemporal patterns even if they are overlapping.
Using simulated data, we tested the performance of the method in various situations. Furthermore, we tested the performance of this method
using real MEG data during a visual stimulation task. We examined
whether the proposed method could estimate the VEFs without using
stimulus onset information. Finally, the method was applied toresting-state fMRI and MEG data. All of the above analyses conﬁrmed
the applicability and usefulness of the proposed method for real brain
activity data.
Proposed method to estimate spatiotemporal patterns
Assumption and purpose
Resting-state data are assumed to contain several unknown spatiotemporal patterns at unknown onsets (Fig. 1A), which is expressed asyðchÞ ðt Þ ¼

Ik
K X
X
k¼1

ðchÞ pk
t−τk;i þ 1 þ vðchÞ ðt Þ;ð1Þ

iwhere y(ch)(t) is resting-state data at channel ch, K is the number ofspatiotemporal patterns, Ik is the number of onsets for the k-th spatio(t) is the k-th spatiotemporal pattern, τk, i is the
temporal pattern, p(ch)ki-th onset of the k-th spatiotemporal pattern, and v(ch)(t) is noise atchannel ch. By introducing onset time seriesuk ðt Þ ¼

1 t ∈ τ k;1 ; ⋯; τ k;Ik
0 Otherwise;

ð2Þ

Onset

ChannelTime

Spatiotemporal
pattern 1

Time
Spatiotemporal
pattern 2
Spatiotemporal
pattern 3
Resting-statedata

B

Spatiotemporal pattern(ch)p1 (t)

Onset time seriesNoise

u1 (t)

v

(ch)

(t)
Resting-state data

N

+

y

(ch)

(t)

(ch)

pK (t)uK (t)

Fig. 1. Assumption of proposed method (STeP). (A): Resting-state data is assumed to contain several unknown spatiotemporal patterns at unknown onsets. Spatiotemporal patterns aredeﬁned as activities represented by two-dimensional matrices of channel × time. (B): Schematic representation of Eq. (3). Note that uk(t) takes a binary (0 or 1) value.

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

u becomes small. Therefore, the objective function is deﬁned as

Eq. (1) can be rewritten in a convolution form asyðchÞ ðt Þ ¼

K X
N
X

253

ðchÞ
pk ðnÞuk ðt−n þ 1Þ þ vðchÞ ðt Þ;ð3ÞRðp; uÞ ¼

"
CH X
T
X

#2
K X
N
X
ðchÞ
yðchÞ ðt Þ−
pk ðnÞuk ðt−n þ 1Þ :ch¼1 t¼1

k¼1 n¼1

ð4Þ

k¼1 n¼1

STeP searches for the optimal p and u that minimize R(p, u).
where N is the length of spatiotemporal patterns (Fig. 1B). Eq. (3) canbe regarded as the ﬁnite impulse response (FIR) model, where p =(t)| k = 1 : K, ch = 1 : CH, t = 1 : N} and u = {uk(t)|k = 1 : K, t = 1 : T}{p(ch)
kcorrespond to the impulse response and the input, respectively. In thefollowing, the notation x = 1 : X is used to represent x = 1 , ⋯ , X. Notethat u takes a binary (0 or 1) value and is unknown in this study. Furthermore, the waveforms of p can be different across channels and arealso unknown.
The purpose of the proposed method is to estimate p and u from y ={y(ch)(t)| ch = 1 : CH, t = 1 : T}, given K and N. Hereafter, this method iscalled STeP (SpatioTemporal Pattern estimation).
Objective function
STeP attempts to estimate p and u so that the power of the residual
error between observed and reconstructed data by the estimated p andAOptimization algorithm
If u is given, p can be estimated by the least squares method. On theother hand, if p is given, we can search for u. Therefore, from an initial u,
we can obtain an estimate of u by alternately iterating the updates of p
and u. However, the converged u may be a local minimum of the objective function [Eq. (4)]. To avoid local minima, we need to repeat thisprocedure from various initial values of u. Furthermore, to generate
the initial u, we need to set the number of onsets, but this is unknown.
To deal with these issues, we developed an optimization algorithm
consisting of three hierarchical levels: top, middle, and bottom(Fig. 2A). The top level is for estimating the necessary number of onsets.
The middle level is for avoiding local minima by generating various
initial values of u. The bottom level is for actually estimating u. These
three levels are described below.

Top level
Start

Middle level
Number of onsets
M=2

Bottom level
Generate
M onsets

Generate
M onsets

Set target onset

Estimate onsets
Estimate E

Estimate 1
M=M+2

Estimate
spatiotemporal
patterns

Changetarget onset

Updatespatiotemporal
patterns

Select best result
Update target onset

No

Residual error
stops decreasing?
No
Converge?

Yes

Yes

End

B

C
u and pTop level

residual error

Magnitude of

Resting-state data

Stop increasing M

M

Selected u

Middle level
Initial uM

Converged u

Bottom level

Fig. 2. Algorithm of proposed method (STeP). (A): Flowchart to estimate spatiotemporal patterns and their onsets. At top level, onsets and spatiotemporal patterns are estimated while
gradually increasing the number of onsets M. At middle level, given M, estimation of onsets is repeated for E times using different initial onsets, and the best onsets are selected. At
bottom level, onsets are estimated from given initial onsets. (B): Schematic representation of top level. As M increases, magnitude of residual error becomes smaller until M reaches a
sufﬁcient value. (C): Information ﬂow across the top, middle, and bottom levels.

254

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

Top level
This level is for estimating the necessary number of onsets. We ﬁrst
set the number of onsets for each spatiotemporal pattern to M. Becausesome of the M onsets generated at the middle level can be removed at
the bottom level (see middle and bottom levels), M deﬁnes the upper
limit of the number of onsets for each spatiotemporal pattern. We set
common M across k. Then, we estimate u (see middle and bottom
levels). Once u is obtained, p can be estimated by the least squares
method as follows. Eq. (3) is rewritten in a matrix form asa spatiotemporal pattern when M = 6, now (M = 8) we generate 8 onsets by adding 3 onsets to the 5 onsets. Such randomness is introduced
only at this level for adding the onsets.
Bottom level
At this level, we actually estimate u from the initial u. Namely, onsetsτ = {τk,i | k = 1 : K, i = 1 : Ik} are estimated by sequentially updating each
onset one by one. Once τ is estimated, u can be determined according
to Eq. (2). Let τk;~ ~i denote the target onset to be updated. We iteratep-step Update spatiotemporal patterns using onsets except for thetarget onset τk;~ ~iu-step Update the target onset τ~k;~i so that the residual error becomes
smaller

Y ¼ UP þ V;
where2

yð1Þ ðT Þ
6 yð1Þ ðT−1Þ
6
Y ¼4⋮
yð1Þ ðNÞ

3⋯yðCHÞ ðT Þ
ðCH Þ
⋯ y
ðT−1Þ 77;
5⋯⋮⋯
yðCHÞ ðNÞ2

⋯ u1 ðT−N þ 1Þ
u1 ðT Þ6 u1 ðT−1Þ ⋯u1 ðT−NÞ
U¼6
4
⋮
⋯
⋮
⋯
u1 ð1Þ
u1 ðNÞ

2

ð1Þ

p1 ð1Þ
6
6 ð1Þ⋮
6 p ðNÞ
6 1
P¼6
6 ð1Þ⋮
6 p ð1Þ
6 K
4⋮
ð1Þ
pK ðNÞ

⋯
⋯
⋯
⋯
⋯
⋯
⋯

2

vð1Þ ðT Þ
6 vð1Þ ðT−1Þ
V ¼6
4
⋮
ð1Þ

v ðN Þ

ðCH Þ

p1

ð1Þ

while changing the target onset with the following order: ~k ¼ 1 : K;~i ¼ 1 : IK . The two steps are described below.
At p-step, onset time series not containing the target onset τ~k;~i are3⋯
uK ðT Þ⋯ uK ðT−N þ 1Þ⋯ uK ðT−1Þ ⋯uK ðT−NÞ 7
7;
5
⋯
⋮
⋯
⋮
⋯
uK ðNÞ
⋯uK ð1Þ

3

7
7
ðNÞ 7
7
7;
⋮
7
ðCH Þ
pK ð1Þ 775⋮
ðCHÞ
pK ðNÞ
⋮u0k ðt Þ ¼

0
uk ðt Þ⋯⋯⋯⋯t ¼ τk;~ ~iOtherwise:0

Using u0 ¼ fu k ðtÞjk ¼ 1 : K; t ¼ 1 : Tg, p ′ = {p′k(ch)(t)| k = 1 : K, ch =1 : CH,t = 1 : N} is obtained by Eq. (5).
At u-step, a residual error is calculated by

ðCHÞ

p1

K

N

k¼1

n¼1

r ðchÞ ðt Þ ¼ yðchÞ ðt Þ  ∑ ∑ p0k ðchÞ ðnÞu0k ðt  n þ 1Þ:~spatiotemBecause u′ does not assume the target onset τ ~k;~i , the k-thporal pattern p0~k ðchÞ ðtÞ is expected to remain in r(ch)(t) around τ ~k;~i .

3Therefore, a candidate time point for the target onset is obtained byvðCHÞ ðT ÞvðCHÞ ðT−1Þ 7
7:
5
⋮
ðCH Þ
v
ðN Þ

T
argmin CH∑ ∑t 0 ∈ t ch¼1 t¼1

t new ¼hi2
ðchÞr ðchÞ ðt Þ  p0 ~ ðt  t 0 þ 1Þ ;kwhere t is the set of time points between the previous and next onsets ofthe target onset ½τk;
~ ~i−1 þ 1; ⋯; τ ~k;~iþ1 −1. The target onset is updated to

The least squares solution of P ishi1
^ ¼ UT UPU T Y:

generated bytnew only if
ð5Þ

CH

T

ch¼1

t¼1

∑ ∑

Finally, we calculate R(p, u).
As M increases, R(p, u) becomes smaller until M reaches a sufﬁcientvalue (Fig. 2B). Therefore, we gradually increase M until R(p, u) stops
decreasing. In this study, we increased M from 2 by 2 and then stopped
increasing it if R(p, u) did not exceed its minimum value in three consecutive loops, or when M reached T/N.
Middle level
At this level, given M, we try to ﬁnd the global minimum solution of
u. Because R(p, u) has local minima, the initial u does not need toconverge to the global minimum. To ﬁnd the global minimum, we repeat the estimation of u (see bottom level) for E times from different
initial u values and select the best u that minimizes R(p, u), where p is
obtained using u by Eq. (5). E was set to 30 in this study. The selected
u is preserved for use in generating the initial u next time (i.e. when
M =M + 2).
When M = 2, the initial u is generated by using only random
numbers. Otherwise, it is generated using the preserved u and random
numbers. Suppose that the current and last M are m and m− 2, respectively. We generated the m initial onsets for each spatiotemporal pattern by adding onsets generated by random numbers to the onsets
estimated when M = m − 2. For example, if 5 onsets are estimated for

h
i2
CH
T
r ðchÞ ðt Þ  p0k~ ðchÞ ðt  t new þ 1Þ b ∑ ∑ r ðchÞ ðt Þ2 :
ch¼1

t¼1

Otherwise, the target onset is removed by assuming the onset is notnecessary within t. Therefore, the number of estimated onsets can besmaller than M and different across k.
The decision of convergence is conducted once in updating all of theonsets. We estimate p using all onsets by Eq. (5) and then calculateR(p, u). We exit the loop if R(p, u) is not smaller than that at the previous
decision.
Different spatiotemporal patterns are automatically assigned to different k in estimating u at the middle and bottom levels. This is because
assigning different spatiotemporal patterns to different k minimizesR(p, u) more than the other approaches, such as assigning the same spatiotemporal pattern to different k.
Fig. 2C illustrates the information ﬂow across the top, middle, andbottom levels. At the top level, we set M and pass it to the middle
level. At the middle level, given M, we generate the initial u and passit to the bottom level. At the bottom level, from the initial u, we estimate
u and return it to the middle level. These steps are repeated for E times
with different initial u, and as a result E-estimated u are obtained. At the
middle level, from the E-estimated u, we select the best u and return itto the top level. Furthermore, we preserve it for generating the initial u
when M =M + 2.

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

Materials and methods

255

noise to the signal time series (Fig. 3C). The standard deviation (SD) of
the noise was either 0.22, 0.40, or 0.71 corresponding to the SNRs of 0,−5, and −10, respectively. The SNR is deﬁned as

Simulation test
The performances of STeP were examined using simulated MEG/EEG
data.
Simulated data were generated as follows. Spatiotemporal patterns
were generated while setting the number of channels to 10 (Fig. 3A,left). The lengths of the spatiotemporal patterns and simulated data
were respectively set to 20 and 6000, corresponding to 0.4 and 120 s
when the sampling rate in an ofﬂine analysis is 50 Hz. We generated
150 onsets for each spatiotemporal pattern using random numbers.
From the onsets, the onset time series were generated according to
Eq. (2) (Fig. 3A, right). By summing the convolutions of the spatiotemporal patterns and their onset time series, we obtained a signal time series (Fig. 3B). Simulated data were generated by adding Gaussian whiteCH

10 log10NðchÞT ∑ch¼1 ∑t¼1 pk ðt Þ2CHN ∑ch¼1

T∑t¼1vðchÞ ðt Þ2

;

N
(ch)
(t)2 is the same across k. By applying STeP towhere ∑CH
ch=1 ∑t=1 pk
the simulated data, we estimated the spatiotemporal patterns and theironsets.
Before quantifying the estimation accuracy, it was necessary tochange the assignment of estimated spatiotemporal patterns to k.
This is because different spatiotemporal patterns are arbitrarily
assigned to different k in the estimation. Furthermore, it was also
necessary to adjust the average of the estimated onsets for each spatiotemporal pattern. This is because the shape of a spatiotemporal
pattern shifts forward or backward depending on the deﬁnition of
its onsets, and these are also arbitrarily determined in theA

ChannelTrue spatiotemporal patterns
2
6
10
2
6
10
2
6
10
2
6
10
2
6
10
10

0.5

1

−0.5
0.5

0
1

−0.5
0.5

0
1

−0.5
0.5

0
1

−0.5
0.5

0
1

−0.5

0

True onset time series

100

20

Time

200

300

Time
ChannelB

Signal
2

2
6
10

−1
100

200

300

Time
ChannelC

Simulated data = Signal + Gaussian white noise
2
6
10

2
−1
100

200

300

Time

D

ChannelEstimated spatiotemporal patterns
2
6
10
2
6
10
2
6
10
2
6
10
2
6
10

10

Time

20

0.5

1

−0.5
0.5

0
1

−0.5
0.5

0
1

−0.5
0.5

0
1

−0.5
0.5

0
1

−0.5

0

Estimated onset time series

100

200

300

Time

Fig. 3. Basic simulation test. (A): True spatiotemporal patterns (left) and their onset time series (right). (B): Signal time series generated by summing convolutions of true spatiotemporal
patterns and their onset time series. (C): Simulated data generated by adding Gaussian white noise to signal time series. (D): Estimated spatiotemporal patterns (left) and their onset timeseries (right).

256

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

estimation. We changed the assignment of the estimated spatiotemporal patterns and adjusted the average of the estimated onsets sothat the difference between the true and estimated spatiotemporal
patterns was minimized. This procedure was also conducted whencomparing spatiotemporal patterns estimated from different dataor with different parameters.
Correlation coefﬁcient of spatiotemporal patterns
To quantify the estimation accuracy of the spatiotemporal patterns,
we calculated the correlation coefﬁcients between the true andestimated spatiotemporal patterns and then averaged them across thespatiotemporal patterns.
Normalized distance from true onsets
To quantify the estimation accuracy of the onsets, we calculated two
metrics: normalized distance from the true onsets and normalized
number of estimated onsets. The normalized distance from the true
onsets represents how close the estimated onsets are to the true ones
compared to the mean inter-onset interval (IOI) of the estimated onsets,
and it was calculated byIk
1 K
1∑ min jτk;i  τ^k; j j;∑
K k¼1 ak  I k i¼1
jwhere Ik is the number of true onsets for the k-th spatiotemporal
pattern, τk , i is the i-th true onset of the k-th spatiotemporal pattern,and τ^k; j is the j-th estimated onset of the k-th spatiotemporal pattern.
ak is a normalization value representing the mean IOI and is calculated by
ak ¼

1

^I
k

∑

^I k  1 i¼2τ^k;i  τ^k;i1 ;where ^I k is the number of estimated onsets for the k-th spatiotemporal
pattern. The normalized distance from the true onsets becomes close to
0 if the estimated onsets are close to the true ones.
Normalized number of estimated onsets
The normalized number of estimated onsets represents how many
onsets are estimated compared to the number of true onsets, and it
was calculated by
1 K ^I k:∑
K k¼1 I k

Human Subject Review Committee. All of them had normal orcorrected-to-normal visual acuity.
The experimental design is shown in Fig. 4. A MEG recording session
consisted of three 5-min consecutive periods: Rest 1, Stimulation, andRest 2. Throughout the session, the subjects were instructed to ﬁxate
on a white cross. During Stimulation, a visual stimulus was presented
100 times with a duration of 0.1 s on the left-hand side of the whitecross. The inter-stimulus intervals were 3–4 s.
The MEG were recorded with a whole-head 400-channel system(210-channel Axial Gradiometer and 190-channel Planar Gradiometer;
PQ1400RM; Yokogawa Electric Co., Japan). The sampling frequency was
1 kHz. An electrooculogram (EOG) value was simultaneously recorded.
In ofﬂine analyses, the MEG data were passed through a low-pass FIR
ﬁlter with a cutoff frequency of 8 Hz, sampled at 50 Hz, and passed
through a high-pass FIR ﬁlter with a cutoff frequency of 1 Hz. Using reference sensor data, environmental noise was removed by time-shift
Principal Component Analysis (PCA) (de Cheveigné and Simon, 2007).
EOG artifacts were removed by generating a multiple linear regression
model to predict eye-movement-related components in the MEG datausing the EOG data, and then the prediction was removed from the
MEG data. Cardiac artifacts and sensor noise were removed by ICA(Jung et al., 2001).
To examine whether STeP could estimate the VEF without using
stimulus onset information, we applied STeP to the continuous MEG
data during Stimulation. We used 10 axial channels where the VEF
was large, and we set the length of spatiotemporal patterns to 0.5 s
based on the waveform of the VEF. For each subject, we determined
the number of spatiotemporal patterns based on reproducibility of estimation results. After estimating spatiotemporal patterns, we selected
the spatiotemporal pattern corresponding to the VEF and adjusted theaverage of its estimated onsets so that the difference between the VEF
and the selected spatiotemporal pattern became smallest. To quantify
the similarity between the VEF and the selected spatiotemporal pattern,
we calculated the correlation coefﬁcient between them. To test the statistical signiﬁcance of the correlation coefﬁcient, we generated 1000
surrogate values for the correlation coefﬁcient using IOI randomized onsets and calculated p = Ns/1000, where Ns is the number of surrogate
values larger than the actual correlation coefﬁcient. To quantify how
many stimulus onsets were detected, we calculated the detection rateby Nd/100, where Nd is the number of trials in which the onsets were
estimated between the stimulus onsets ± 0.02 s. The statistical signiﬁcance of the detection rate was tested in the same way as the correlation
coefﬁcient. Furthermore, to quantify the extent to which the VEF isdominant in the MEG data, we calculated the contribution ratio of the
VEF byThis value becomes larger than 1 if there are false positive onsets.
10

MEG

1

To show the applicability of STeP for real MEG data, we conducted a
MEG experiment.
Eleven healthy subjects (ages 26.8 ± 7.5 years [mean ± SD]) participated in this experiment. All subjects gave written informed consent
for the experimental procedures, which were approved by the ATR

T∑ch¼1 ∑t¼1 yðchÞ ðt Þ2
10
∑ch¼1

T

∑t¼1 yðchÞ ðt Þ2;

where y -(ch)(t) was obtained by subtracting the VEF from the MEG data.
To estimate repetitive spatiotemporal patterns in resting-state MEG
data, we applied STeP to the MEG data during Rest 1. The MEG data

Rest 1

Stimulation

Rest 2(5 min)(5 min)(5 min)

+

Beep

+

+

+

0.1 s

3–4 s

0.1 s

Beep

+

Time
Fig. 4. Experimental design of MEG experiment. A MEG recording session consists of three consecutive periods: Rest 1, Stimulation, and Rest 2. Throughout the experiment, subjects were
instructed to ﬁxate on a white cross.

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

were preprocessed as described above except that the cutoff frequencyof the low-pass ﬁlter was 25 Hz. We used 210 axial channels.
fMRI
To show the applicability of STeP for fMRI data, we conducted a
resting-state fMRI experiment.
Ten healthy subjects (ages 29.2 ± 7.6 years) participated inthis experiment. All gave written informed consent for the experimental procedures, which were approved by the ATR Human
Subject Review Committee, and all had normal or corrected-tonormal visual acuity.
The experiment was done in a 10-min resting-state condition. The
subjects were instructed to ﬁxate on a white dot, to let their mind
wander, and to not focus on any one thing.
Three Tesla MR scanner (MAGNETOM Trio 3 T; Siemens, Germany)
was used to obtain the structural and functional MRI data. The following
are the acquisition parameters for the T1-weighted images: repetition
time 2300 ms, time of echo 2.98 ms, ﬂip angle 9°, slice thickness
1 mm, ﬁeld of view 256 × 256 mm, and imaging matrix 256 × 256
with 240 slices. The following are the acquisition parameters for the
echo-planar images (EPIs): repetition time 2500 ms, time of echo
30 ms, ﬂip angle 80°, slice thickness 4 mm, ﬁeld of view
211.84 ×211.84 mm, and imaging matrix 64 ×64 with 40 slices.
In ofﬂine analyses, the fMRI data were preprocessed by SPM8
(Welcome Department of Cognitive Neurology, UK). Here, head motion
and slice-timing were corrected. The images were spatially normalized
to match the MNI template and smoothed with an 8-mm full-width at
half-maximum (FWHM) Gaussian ﬁlter. The time series of each voxel
was high-pass ﬁltered to 1/128 Hz. To remove motion artifacts, we
subtracted components correlated with head motion from the timeseries of each voxel. After normalizing the time series of each voxel toobtain a mean of 0 and an SD of 1, the preprocessed fMRI data of thegray matter for all subjects were concatenated together. To verify the effect of preprocessing, we applied both ICA and seed-based correlation
analysis to the concatenated data. This revealed spatial patterns
consistent with the DMN (Fox et al., 2005; Smith et al., 2009; Tonget al., 2015) (Supplementary material), suggesting the effectiveness ofthe preprocessing.

B

A
True
2

0.5

10
2

−0.5
0.5

Channel

−0.5
0.5

0.1

10
2

0.4

−0.5

10
2

0.5

0.4

−0.5

0

10
2

0.5

10
2

−0.5

10

0.4
0.2
0.4

6

6
10
10

Time

20

D
0

2

0
10

Time

20

Olshausen‘smethod (2003)
0.2

2
6

−0.1

10
2

0.1
0

10
2

0

6

6

6
10
2

The performance of STeP was tested using simulated data.
We generated ﬁve spatiotemporal patterns (Fig. 3A, left) and generated their onsets using random numbers (Fig. 3A, right). By summing
the convolutions of the spatiotemporal patterns and their onset time series, we obtained a signal time series (Fig. 3B). We can see that the spatiotemporal patterns are overlapping (Fig. 3B). Simulated data weregenerated by adding Gaussian white noise to the signal time series(Fig. 3C). The SD of the noise was set to 0.40 so that the SNR would become −5. The simulated data are too noisy to visually distinguish thespatiotemporal patterns and their onsets (Fig. 3C). We estimated thespatiotemporal patterns and their onsets by applying STeP to the simulated data (Fig. 3C). In the estimation, we set the number and length of
the spatiotemporal patterns to the true values of 5 and 20, respectively.
The estimation required about 8 min on a Xeon processor (3.2 GHz × 8
cores).
Fig. 3D shows the estimated spatiotemporal patterns and onsets. The
estimated spatiotemporal patterns resemble the true ones (Fig. 3A andD, left) (the correlation coefﬁcient between the true and estimated spatiotemporal patterns was 0.98). The estimated onsets are close to thetrue onsets (Fig. 3A and D, right), with the normalized distance fromthe true onsets at 0.02. The number of estimated onsets is almost thesame as that of true onsets (Fig. 3A and D, right), with the normalized
number of estimated onsets at 1.18. These results indicate that STeP successfully estimated the spatiotemporal patterns and their onsets fromthe noisy simulated data (Fig. 3C), even though the spatiotemporal patterns were overlapping.

6
0

10
2
6

6

Basic simulation test

6

6

6
10
2

Results

ICA
0.5

2

We then applied STeP to the concatenated fMRI data. To shorten thecomputation time, we used PCA to reduce the dimension of the fMRI
data from the number of voxels to that of the samples while keeping
its amplitude information. From this low-dimensional data, we estimated the onsets of spatiotemporal patterns. Then, using the estimated onsets, we estimated spatiotemporal patterns from the original fMRI data
by Eq. (5).

C

Majeed et al.‘smethod (2011)

6

6

257

0

10
2

0.2

6

0

10
2

0.1
0
−0.1
0.2

6
−0.1

10
2

0.1

10
2

6

0

6

0

10
2

−0.1
0.05

10
2

0.2

6

0

−0.1

10

6
10
10

Time

20

10

20

Time

Fig. 5. Performances of three other methods. (A): True spatiotemporal patterns, which are the same as those shown in Fig. 3A, left. (B): Majeed et al.'s method (2011). Five representative
templates are shown. (C): ICA. Mixing vectors of ICs were reorganized to form spatiotemporal patterns. (D): Olshausen's method (2003). Estimated basis functions are shown.

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

Testing other methods
To clarify the differences between STeP and other approaches, three
existing methods were also applied to the same simulated data
(Fig. 3C). Their performances were also quantiﬁed by the correlation coefﬁcients between the true spatiotemporal patterns (Fig. 5A) and the estimated ones, which were ordered and shifted along time to maximize
the correlation coefﬁcients.
First, we tested Majeed et al.'s method (2011), which is basically atemplate-matching algorithm. Fig. 5B shows examples of the estimatedtemplates. The estimated templates do not resemble the true spatiotemporal patterns (Fig. 5A and B) (the correlation coefﬁcient between
the true spatiotemporal patterns and the templates was 0.50). This
performance is attributable to the fact that, in the simulated data, thespatiotemporal patterns are overlapping (Fig. 3B), and this method
cannot separate them. In the case of STeP, overlapping spatiotemporal
patterns can be separated (Fig. 3D, left).
Second, we tested ICA. To obtain spatiotemporal patterns by ICA, we
ﬁrst generated time-shifted data. We increased the dimension of thesimulated data from CH to CH × N by treating the data at time t − n(n = 1 , ⋯ , N − 1) as additional channels. Then, we applied ICA to thetime-shifted data while regarding the time-shifted data at each time
point as a sample and setting the number of ICs to 5. We used runica.m
from EEGLAB, version 13.4.4b (Delorme and Makeig, 2004). The
resulting mixing vectors of the ICs were reorganized to form spatiotemporal patterns ([CH,N]). Fig. 5C shows the estimated mixing vectors. The
estimated mixing vectors do not resemble the true spatiotemporal
patterns (Fig. 5A and C) (the correlation coefﬁcient between the true
spatiotemporal patterns and the estimated mixing vectors was 0.53).
This is attributed to the fact that ICA is not designed for estimating repetitive spatiotemporal patterns.
Finally, we tested Olshausen's method (2003), which has a similar
generative model and purpose as STeP. Its main difference from STeP
is that u [Eq. (3)] is assumed to represent coefﬁcients of basis functions
p and takes continuous values. We applied his method to the simulated
data (Fig. 3C) while setting the number of basis functions to 5. Fig. 5Dshows the estimated basis functions. The estimated basis functions do
not resemble the true spatiotemporal patterns (Fig. 5A and D) (thecorrelation coefﬁcient between the true spatiotemporal patterns andthe estimated basis functions was 0.54). This is attributable to the factthat the estimated u took continuous values as described above (notshown), although the true u took binary values (Fig. 3A, right). In thecase of STeP, u is assumed to represent the onset timing of spatiotemporal patterns and should be a binary (0 or 1) value. This works as a strong
constraint and improves the estimation accuracy of STeP.
Detailed simulation tests
We examined the performance of STeP in more detail. We repeatedly evaluated the estimation accuracies of spatiotemporal patterns andtheir onsets at the SNRs of − 10, − 5, and 0. The simulated data were

1

−10

−5

SNR

0

0.2

0.1

0
−10

−5

SNR

0

2
Accurate

0.6

C

Accurate

0.8

Normalized distancefrom true onsets

B

Accurate

Correlation coefficient of
spatiotemporal patternsA

generated with the same parameters as the ﬁrst simulation test
(Fig. 3) except for the SNRs. In the estimation, we set the number andlength of the spatiotemporal patterns to the true values of 5 and 20, respectively. Simulation tests were conducted in 20 runs using different
onsets and noise.
Fig. 6 shows the estimation accuracies quantiﬁed by the three
metrics: the correlation coefﬁcients of the spatiotemporal patterns(A), the normalized distances from the true onsets (B), and the normalized numbers of estimated onsets (C). When the SNR is 0, all three metrics show high estimation accuracies for all runs. This result indicates
that, when the SNR is sufﬁciently high, STeP robustly estimates the spatiotemporal patterns and their onsets accurately. When the SNR is −10,
the normalized numbers of estimated onsets are much larger than 1.
This result indicates that, when the SNR is low, STeP estimates many
false positive onsets.
In actual application, the number of spatiotemporal patterns is unknown, but we need to set it. Here, we examined the performance of
STeP when the assumed number of spatiotemporal patterns waswrong. The simulated data were generated with the same parameters
as the ﬁrst simulation test (Fig. 3). In the estimation, we set the numberof spatiotemporal patterns to either 3, 4, 5 (true), 6, or 7. The simulation
tests were conducted in 2 runs using different onsets and noise. Fig. 7Ashows the estimated spatiotemporal patterns in the ﬁrst run. When theassumed number of spatiotemporal patterns is smaller than true (b 5),
some of the estimated spatiotemporal patterns seem to be contaminated with other spatiotemporal patterns. When the assumed number ofspatiotemporal patterns is larger than true (N 5), some of the estimatedspatiotemporal patterns seem to have low SNR (Fig. 7A). To examine
the reproducibility of the estimated spatiotemporal patterns, we calculated the correlation coefﬁcients of the estimated spatiotemporal patterns between the ﬁrst and second runs. Fig. 7B shows the correlation
coefﬁcients. When the assumed number of spatiotemporal patterns islarger than true (N 5), some of the estimated spatiotemporal patterns
show low correlation coefﬁcients (b0.6), indicating that they had low
reproducibility. This result suggests a way to determine the number ofspatiotemporal patterns: examine the reproducibility of estimated spatiotemporal patterns by dividing data into two parts.
The length of spatiotemporal patterns is also unknown, but we needto set it. We examined the performance of STeP when the assumed
length of spatiotemporal patterns was wrong. The simulated data
were generated with the same parameters as the ﬁrst simulation test
(Fig. 3). In the estimation, we set the length of the spatiotemporal patterns to either 10, 15, 20 (true), 25, 30, 35, or 40. The simulation tests
were conducted in 20 runs using different onsets and noise. Fig. 8Ashows the correlation coefﬁcients between the true and estimated spatiotemporal patterns. The correlation coefﬁcients at the lengths of 15,
25, and 30 are not different from those at the length of 20 (true)(p N 0.05, two-tailed Wilcoxon signed-rank test, Bonferroni-corrected).
This result indicates that the estimation accuracy did not greatly decrease unless the assumed length of the spatiotemporal pattern was

Normalized number of
estimated onsets

258

1.5

1

−10

−5

0

SNR

Fig. 6. Detailed performance of STeP. (A): Correlation coefﬁcients between true and estimated spatiotemporal patterns. (B): Normalized distances from true onsets. (C): Normalized
numbers of estimated onsets. In all ﬁgures, error bars represent SDs across runs.

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

A

B

Estimated spatiotemporal patterns

259

Reproducibility of
estimated spatiotemporal patterns
1
0.8

2
3

0.6

4
0.4

5
6

True

0.2

True

0.5

7

−0.5
3

4

5

6

3

7

4

5

6

Correlation coefficient

Spatiotemporal pattern

1

0

7

Assumed number of spatiotemporal patterns
Fig. 7. Performance of STeP when assumed number of spatiotemporal patterns is wrong. (A): Estimated spatiotemporal patterns for each assumed number of spatiotemporal patterns. True
number of spatiotemporal patterns is 5. (B): Reproducibility of estimated spatiotemporal patterns. Correlation coefﬁcients of estimated spatiotemporal patterns between ﬁrst and second
runs are shown.

far from the true value. Therefore, we need not know the exact length of
the spatiotemporal patterns.
In the case of real brain activity data, it is possible that the amplitudes of spatiotemporal patterns vary from onset to onset. To examine
the performance of STeP for such data, we applied STeP to simulateddata, including spatiotemporal patterns whose amplitudes are variable.
The simulated data were generated with the same parameters as theﬁrst simulation test (Fig. 3), except that the spatiotemporal patterns
were multiplied by 1 +w, where w is a Gaussian white noise that is different across the onsets and spatiotemporal patterns. The SDs of w were
set to 0, 0.2, 0.4, 0.6, 0.8, and 1. In the estimation, we set the number andlength of the spatiotemporal patterns to the true values of 5 and 20, respectively. The simulation tests were conducted in 20 runs for each SDof w using different onsets and noise. Fig. 8B shows the correlation coefﬁcients between the true and estimated spatiotemporal patterns as afunction of the SDs of w. The correlation coefﬁcients at SD = 0.2, 0.4,
0.6, 0.8, and 1 are not different from those at SD = 0 (p N 0.05, twotailed Wilcoxon signed-rank test, Bonferroni-corrected). This result indicates that the variability of the amplitudes did not greatly decrease
the estimation accuracy. This is because, as the variability of the amplitudes becomes large, huge spatiotemporal patterns sometimes appear
and STeP easily detects them while missing small ones. For the VEF
shown in Fig. 9C, the SD of w was estimated at 0.36. Considering these

Effect of assuming wrong length of
spatiotemporal patterns
1

0.9

True

0.8
10

20

Performance test with real MEG data
We examined the performance of STeP with real MEG data. It is generally believed that repeated presentation of a visual stimulus evokes astereotypical MEG pattern, the so-called VEF. STeP is proposed to estimate such stereotypical spatiotemporal patterns without using theironset information. Therefore, we tested whether STeP could estimate
the VEF without using the stimulus onset information.
We applied STeP to the continuous MEG data during Stimulation(Fig. 4) at 10 axial channels where the VEFs are large. We set the lengthof spatiotemporal patterns to 0.5 s based on the waveforms of the VEFs.
For each subject, we determined the number of spatiotemporal patterns
based on the reproducibility of the estimation results. Fig. 9 shows theresults of one subject. Fig. 9A shows the 10 selected channels. To determine the number of spatiotemporal patterns, we divided the MEG data
into two parts, each of which was 2.5-min data, and estimated spatiotemporal patterns separately while setting the number of spatiotemporal patterns to either 1, 2, 3, 4, or 5. Fig. 9B shows the correlation
coefﬁcients between the spatiotemporal patterns separately estimatedfrom the two parts. When the assumed number of spatiotemporal patterns is larger than 1, some of the estimated spatiotemporal patterns

B
Correlation coefficient

Correlation coefficientA

results, the variability of the amplitudes of spatiotemporal patterns is
considered no problem for STeP.

30

40

Assumed length of spatiotemporal patterns

Effect of variable amplitude of
spatiotemporal patterns
1

0.9

0.8
0

0.2

0.4

0.6

SD of amplitude

0.8

1
p<0.05

Fig. 8. Performance of STeP when our assumption is not correct. (A): Performance of STeP when assumed length of spatiotemporal patterns is wrong. True length of spatiotemporal
patterns is 20. Correlation coefﬁcients between true and estimated spatiotemporal patterns are shown. (B): Performance of STeP when amplitudes of spatiotemporal patterns arevariable. Correlation coefﬁcients between true and estimated spatiotemporal patterns are shown. Horizontal axis represents SD of w, where spatiotemporal patterns were multipliedby 1+w. In both ﬁgures, error bars represent SDs of correlation coefﬁcients across runs.

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

A

B

Reproducibility of estimatedspatiotemporal patterns

Spatiotemporal
pattern260

1

Selected channels

1

0.6

3

0.4

4

0.2

5
1

2

3

4

5

Correlation
coefficient

0.8

2

0

Assumed number of spatiotemporal patterns

D

VEF

x 10

Channel

Trial

2
0

Channel

x 10
2
1
0

0.2

−1
[T]

0.4

40
80

−2
[T]

Spatiotemporal pattern

Estimated onsets

−13

−13

Frequency

C

60
40
20
−2

−1

Time [s]012Time [s]E

MEG data

Channel

x 10
4

−13

2
0
121

122

123

−2
[T]

Time [s]Stimulus onset
Period of estimated spatiotemporal patternFig. 9. Estimation accuracy of a subject's VEF. (A): Selected 10 channels. (B): Reproducibility of estimated spatiotemporal patterns. Correlation coefﬁcients between spatiotemporal
patterns separately estimated from the two parts of the divided data are shown. Based on this result, we set the number of spatiotemporal patterns to 1 for this subject. (C): VEF (top)and estimated spatiotemporal pattern (bottom). (D): Raster plot (top) and histogram (bottom) of estimated onsets. Time 0 corresponds to stimulus onsets. (E): Example of MEG data
at 10 selected channels. Red triangle indicates a stimulus onset. Black bars indicate periods of estimated spatiotemporal pattern, onsets of which correspond to estimated onsets. In(C) and (E), locations of channels are indicated by color in (A).

show low correlation coefﬁcients (b 0.7), that is, these patterns have low
reproducibility. Therefore, we set the number of spatiotemporal patterns to 1 for this subject. Fig. 9C shows the VEF (top), which was obtained by averaging the MEG data triggered at the stimulus onsets,
and the estimated spatiotemporal pattern (bottom) obtained by setting
the number of spatiotemporal patterns to 1. The estimated spatiotemporal pattern resembles the VEF (the correlation coefﬁcient between
them is 0.90, p b 0.01). Fig. 9D shows the raster plot (top) and histogram(bottom) of the estimated onsets, which were generated by segmenting
the estimated onset time series triggered at the stimulus onsets. A large
number of onsets were estimated around time 0, corresponding to thestimulus onsets. The detection rate of the stimulus onsets is 0.92
(p b 0.01), indicating that STeP detected the stimulus onsets at a signiﬁcantly high rate. These results indicate that STeP successfully estimatedthe VEF without using the stimulus onset information.
Actually, many onsets were estimated far from the stimulus onsets(Fig. 9D). In fact, at the estimated onsets far from the stimulus onsets,
there seems to be a spatiotemporal pattern resembling the estimatedone (Fig. 9E, see MEG data under black bars). Han et al. (2008) showedthat visually evoked cortical activity reverberates in subsequent spontaneous waves. The estimated onsets far from the stimulus may reﬂect thereverberation.
The estimation results of all subjects are summarized in Fig. 10.
These ﬁgures show scatter plots of the contribution ratios of the VEFs
and the two metrics, i.e. the correlation coefﬁcients between the VEFs
and the estimated spatiotemporal pattern (A), as well as the detection
rates of the stimulus onsets (B). Each dot corresponds to each subject,and its color represents the number of spatiotemporal patterns determined based on the reproducibility of the estimation results. The correlation coefﬁcients and the detection rates are higher than 0.8 for the
subjects whose contribution ratios of the VEFs are larger than 0.1. This
indicates that STeP successfully estimated the VEFs and their onsets, unless the VEFs were too small. It is assumed that, when the VEFs were too
small, STeP estimated more dominant spatiotemporal patterns than the
VEFs.
In summary, STeP can estimate the VEF without using any stimulus
onset information, indicating the applicability of STeP for real MEG data.

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

B
1

1

0.8

0.8

Detection rate

Correlation coefficient

A

261

0.6
0.4
0.2
0

Number of
spatiotemporal
patterns
1
2
3
4

0.6
0.4
0.2

0

0.1

0.2

0.3

0

0.4

p<0.01
0

0.1

0.2

0.3

0.4

Contribution ratio of VEF
Fig. 10. Estimation accuracies of all subjects' VEFs. (A): Scatter plots of contribution ratios of VEFs and correlation coefﬁcients between VEFs and estimated spatiotemporal patterns.
(B): Scatter plots of contribution ratios of VEFs and detection rates of stimulus onsets. In both ﬁgures, each dot corresponds to one subject, and its color represents the number ofspatiotemporal patterns determined based on the reproducibility of estimation results. Signiﬁcantly large values are indicated by asterisks.

Application to resting-state fMRI data
To demonstrate the usefulness of STeP, we applied it to the restingstate fMRI data.
To determine the number of spatiotemporal patterns, we ﬁrst examine the reproducibility of estimated spatiotemporal patterns. We divided the resting-state fMRI data into two parts, each of which contained
ﬁve subjects' data, and estimated spatiotemporal patterns separately
while setting the number of spatiotemporal patterns to either 1, 2, 3,
4, or 5. The length of spatiotemporal patterns was set to 10 s. Fig. 11A
shows the correlation coefﬁcients between the spatiotemporal patterns
separately estimated from the two parts. When the assumed number ofspatiotemporal patterns is larger than 2, some of the estimated spatiotemporal patterns show low correlation coefﬁcients (b 0.7), that is,these patterns have low reproducibility. Therefore, we decided to set
the number of the spatiotemporal patterns to 2. Furthermore, we
assessed the sensitivity of the estimation results to the assumed lengthof the spatiotemporal patterns. From the whole resting-state fMRI data,
which contains ten subjects' data, we estimated the spatiotemporal patterns while setting the length of the spatiotemporal patterns to either 5,
7.5, 10, 12.5, or 15 s. The number of spatiotemporal patterns was set to
2. Fig. 11B shows the pairwise correlation coefﬁcients between the estimated spatiotemporal patterns of different lengths. When the length of
the spatiotemporal pattern is longer than 5 s, the correlation coefﬁcientsare constantly high (N 0.8). This indicates that the estimation results

Reproducibility of
estimated spatiotemporal patterns

A

were not sensitive to the assumed length of the spatiotemporal pattern
if it was larger than 5 s.
Fig. 12A shows the spatiotemporal patterns obtained by setting thenumber and length of the spatiotemporal patterns to 2 and 10 s, respectively. The two spatiotemporal patterns show positive and negative activities, respectively. In both patterns, the spatial patterns at 10 s seemto reﬂect large activities at veins, suggesting that these spatial patterns
reﬂect the cerebral venous drainage. Furthermore, in both patterns, regions belonging to the same resting-state networks (RSNs) show large
activities at the same time points: 5 s for auditory, 7.5 s for visual andcerebellum (Smith et al., 2009). This is reasonable because an RSN consists of co-activating regions. Tong et al. (2012, 2015) showed that systemic low-frequency oscillations (sLFOs) travel through the entire brain
and reach different voxels with different delays. The propagation pattern in the estimated spatiotemporal patterns (Fig. 12A), from the auditory cortex to veins via the visual and cerebellum cortices, seems to beconsistent with the delay map of sLFOs (Fig. 3A in Tong et al., 2015).
For comparison, we also applied Majeed et al.'s method (2011) tothe same resting-state fMRI data. The length of the template was set
to 10 s. Fig. 12B shows the template estimated by Majeed et al.'s
method (2011). The spatial pattern at 10 s in the template resembles
the spatial pattern at 5 s in the ﬁrst spatiotemporal pattern (the correlation coefﬁcient between them was 0.87) (Fig. 12A, left and B). The spatial pattern at 2.5 s in the template resembles the spatial pattern at 7.5 s
in the second spatiotemporal pattern (the correlation coefﬁcient

Consistency of
estimated spatiotemporal patterns

B

2
0.6
3
0.4
4
0.2
5
1

2

3

4

5

0

Assumed number of spatiotemporal patterns

0.45

5

0.89

0.55

0.62

0.98

0.92

0.87

0.84

0.93

0.8
0.6
10

0.4
0.94
0.2

Correlation coefficient

0.8

Assumed length of
spatiotemporal patterns [s]1Correlation coefficient

Spatiotemporal pattern

1
1

15
5

10

15

0

Assumed length of spatiotemporal patterns [s]Fig. 11. Effects of assumed number and length of spatiotemporal patterns on estimation results from resting-state fMRI. (A): Effects of assumed number of spatiotemporal patterns on
reproducibility of estimation results. Correlation coefﬁcients between spatiotemporal patterns separately estimated from the two parts of the divided data are shown. (B): Effects ofassumed length of spatiotemporal patterns on estimation results. Pairwise correlation coefﬁcients between estimated spatiotemporal patterns of different lengths are shown.

262

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

A

B
Template estimated by
Majeed et al.’s method (2011)

Estimated spatiotemporal patterns
5s

10 s

5s

10 s

5s

10 s

Superior

Inferior

-1.07

1.07

-0.89

0.89

-0.87

0.87

Fig. 12. Spatiotemporal patterns estimated from resting-state fMRI data concatenated across all subjects. (A): Two spatiotemporal patterns estimated by STeP. (B): Template estimated by
Majeed et al.'s method (2011). In these ﬁgures, activities whose absolute values are over 0.3 of maximum value across voxels are shown.

between them was 0.83) (Fig. 12A, right and B). These similarities suggest the validity of the estimated spatiotemporal patterns.
Application to resting-state MEG data
We also applied STeP to the resting-state MEG data. Here, we show
the estimation results of a subject.
To determine the number of spatiotemporal patterns, we ﬁrst examine the reproducibility of estimated spatiotemporal patterns. We divided the resting-state MEG data into two parts, each of which was 2.5-min
data, and estimated spatiotemporal patterns separately while setting
the number of spatiotemporal patterns to either 1, 2, 3, 4, or 5. The
length of spatiotemporal patterns was set to 0.5 s. Fig. 13A shows thecorrelation coefﬁcients between the spatiotemporal patterns separately
estimated from the two parts. When the assumed number ofReproducibility of
estimated spatiotemporal patternsA

spatiotemporal patterns is larger than 1, all of the estimated spatiotemporal patterns show low correlation coefﬁcients (b 0.7). Therefore, we
decided to set the number of spatiotemporal patterns to 1. Furthermore,
we assessed the sensitivity of the estimation results to the assumed
length of the spatiotemporal pattern. From the whole resting-state
MEG data, we estimated the spatiotemporal patterns while setting thelength of the spatiotemporal patterns to either 0.2, 0.3, 0.4, 0.5, 0.6,
0.7, or 0.8 s. The number of spatiotemporal patterns was set to 1. Fig.
13B shows the pairwise correlation coefﬁcients between the estimatedspatiotemporal patterns of the different lengths. When the length of thespatiotemporal pattern is longer than 0.2 s, the correlation coefﬁcientsare constantly high (N 0.7). This indicates that the estimation results
were not so sensitive to the assumed length of the spatiotemporal pattern if it was larger than 0.2 s.

Consistency of
estimated spatiotemporal patterns

B1

Spatiotemporal pattern

3
0.4
4
0.2
5
1

2

3

4

5

0

Assumed number of spatiotemporal patterns

0.71 0.71 0.65 0.63 0.60 0.55

0.2

0.92 0.81 0.92 0.75 0.80
0.84 0.94 0.80 0.83

0.4

0.8
0.6

0.86 0.91 0.92
0.74 0.92

0.6

0.86

0.4
0.2

Correlation coefficient

0.6

Correlation coefficient

0.8
2Assumed length of
spatiotemporal patterns [s]

1
1

0.8
0.2

0.4

0.6

0.8

0

Assumed length of spatiotemporal patterns [s]Fig. 13. Effects of assumed number and length of spatiotemporal patterns on estimation results from resting-state MEG. (A): Effects of assumed number of spatiotemporal patterns on
reproducibility of estimation results. Correlation coefﬁcients between spatiotemporal patterns separately estimated from the two parts of the divided data are shown. (B): Effects ofassumed length of spatiotemporal patterns on estimation results. Pairwise correlation coefﬁcients between spatiotemporal patterns of different lengths are shown.

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

A

263

0.02 s0.04 s

0.06 s0.08 s0.1 s0.12 s0.14 s0.16 s0.18 s0.2 s0.22 s0.24 s0.26 s0.28 s0.3 s0.32 s0.34 s0.36 s0.38 s0.4 s0.42 s0.44 s

0.46 s0.48 s

0.5 s

x 10

−13

1
−1 [T]

B
Cluster 1

x 10

Cluster 2

−14

10

x 10
1

−13

Cluster 3

0

0
[T]

x 10
10

−14

0
[T]

−1 [T]

coefficient

CorrelationC
1

Cluster 1
Cluster 2
Cluster 3

0
−1

0.1

0.2

0.3

0.4

0.5

Time [s]Fig. 14. Spatiotemporal pattern estimated from a subject's resting-state MEG data. (A): Spatiotemporal pattern. (B): Centers of clusters estimated by k-means clustering. (C): Correlation
coefﬁcients between spatiotemporal pattern at each time (A) and centers of clusters (B).

Fig. 14A shows the spatiotemporal pattern obtained by setting thenumber and length of the spatiotemporal patterns to 1 and 0.5 s, respectively. The spatiotemporal pattern contains spatial patterns similar to
the movement-related magnetic ﬁelds at 0.20, 0.24, 0.34 and 0.44 s.
For comparison, we also applied k-means clustering (MacQueen,1967) to the same resting-state MEG data. Here, the resting-state MEG
data at each time point was regarded as a sample, and the number ofclusters was set to 3. Fig. 14B shows the centers of the estimated clusters(Clusters 1–3). Fig. 14C shows the correlation coefﬁcients between thespatiotemporal pattern at each time and the centers of the clusters. Allof the clusters have correlation coefﬁcients close to 1 at different time
points (0.44 s for Cluster 1, 0.12 s for Cluster 2, and 0.48 s for Cluster 3),
indicating the consistency between the results of STeP and k-means clustering. On the other hand, the correlation coefﬁcients are not constant butchange with time. This indicates that the spatiotemporal pattern represents the consecutive brain activities that dynamically change withtime, while each cluster represents a snapshot of them.

Discussion
In this study, we proposed the STeP method to estimate repetitive
spatiotemporal patterns from resting-state brain activity data. The performance tests with the simulated data show that STeP can estimate
spatiotemporal patterns accurately without using their onsets, even if
they are overlapping (Fig. 3). Furthermore, our results also show that
STeP works robustly in various situations (Figs. 6 and 8). The performance tests with real MEG data show that STeP can estimate the VEFs
without using stimulus onset information (Figs. 9 and 10). Finally,
STeP was applied to the resting-state fMRI and MEG data. This revealed
informative spatiotemporal patterns, showing how the brain activities
dynamically change with time (Figs. 12 and 14). These results indicate
the applicability and usefulness of STeP for resting-state brain activity
data.

264

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

Methodological considerations
STeP searches for spatiotemporal patterns and their onsets that minimize the objective function [Eq. (4)]. This is a difﬁcult optimization
problem because there are many local minima. To solve this optimization problem, we developed the algorithm shown in Fig. 2. In this algorithm, each onset is updated sequentially (Fig. 2A, bottom level), so it istime-consuming. However, it is also possible to update all of the onsets
together. For example, deconvolving spatiotemporal patterns out of
resting-state data provides the information of all onsets at once. Actually, in our simulation tests (not shown) this deconvolving algorithm wasvery fast. However, its estimation accuracy was low because it was always trapped in local minima. Other algorithms that update all of theonsets together have shown similar results. Considering both thecomputational cost and the optimization performance, the proposed
algorithm (Fig. 2) is empirically the best among all those we tested.
The performance tests with the simulated data show that, when SNR
is low, estimated onsets include many false positives (Fig. 6C); in otherwords, overﬁtting occurs. A solution to this problem could be adding a
regularization term to the objective function [Eq. (4)]. Imposing a penalty for increasing the number of onsets by a regularization term may
suppress the false positive onsets and increase the estimation accuracy.
Overﬁtting also occurs when the assumed number of spatiotemporal
patterns is larger than the actual value (Fig. 7A). In this study, we suggest a solution to this problem: determine the number of spatiotemporal patterns based on the reproducibility of estimated spatiotemporal
patterns (Fig. 7B). On the other hand, if our model [Eq. (3)] is
reformulated with a probabilistic model, we could avoid overﬁtting
using well-known criteria, such as Akaike's information criterion (AIC)(Akaike, 1974), the Bayesian information criterion (BIC) (Schwarz,1978), and the free energy.
Reliability of estimation results
STeP estimates repetitive spatiotemporal patterns and their onsets ifthey exist. However, even if they do not exist, STeP is forced to extract
spatiotemporal patterns so that the residual error [Eq. (4)] becomes
small. Therefore, it is necessary to test the reliability of estimated spatiotemporal patterns and onsets.
The reliability of estimated spatiotemporal patterns could be tested
by examining their reproducibility. As shown in this study (Figs. 7B,
9B, 11A, and 13A), reproducibility can be examined by estimating the
spatiotemporal patterns from two separate sets of data and assessing
the similarity between them. If reproducibility is high, the estimatedspatiotemporal patterns could be considered reliable. Furthermore,
the reproducibility of estimated spatiotemporal patterns could be acriterion for determining the number of spatiotemporal patterns(Figs. 7B, 9B, 11A, and 13A).
The reliability of estimated onsets could be tested by examining thedata at estimated onsets (see Fig. 9E). If the data at an estimated onset issimilar to the corresponding spatiotemporal pattern, the estimatedonset could be considered reliable. This similarity can be quantiﬁed by
the correlation coefﬁcients.
Relation to other methods
For spike and fMRI data, template-matching algorithms have been
applied to estimate repetitive spatiotemporal patterns (Ikegaya et al.,
2004; Majeed et al., 2011). The template-matching algorithms use asegment starting at a random time point as a template. Therefore, thesegment must contain a pure spatiotemporal pattern that is not overlapping with others. Otherwise, the template-matching algorithms estimate spatiotemporal patterns contaminated by each other (Fig. 5B).
Because spatiotemporal patterns are assumed to overlap in the case of
MEG/EEG data, the template-matching algorithms are not suitable forthese data.

ICA (Beckmann et al., 2005) and clustering methods, such as kmeans clustering (MacQueen, 1967), are very powerful tools to analyze
resting-state brain activity data. However, these methods do not provide spatiotemporal patterns, which are two-dimensional matrices ofchannel × time (see Fig. 1A). Therefore, these methods do not meet
our needs.
If ICA or k-means clustering is applied to time-shifted data, which
are obtained by treating the data at time t − n (n = 1 , ⋯ , N − 1) asadditional channels, the resulting mixing matrices of ICs or the centersof clusters can be regarded as spatiotemporal patterns. However,mixing matrices of ICs or centers of clusters do not necessarily reﬂect
actual spatiotemporal patterns (Fig. 5C). This is because neither ICA
nor k-means is designed to estimate repetitive spatiotemporal patterns.
Therefore, applying ICA or k-means clustering to time-shifted data is notsuitable for our purpose.
The methods of Olshausen (2003) and Anemüller et al. (2003) have
similar generative models and purposes as ours. STeP is different fromthese methods mainly in the assumption of u [Eq. (3)]. In thesemethods, u takes continuous values, assuming that u represents coefﬁcients of bases. In STeP, on the contrary, u takes a binary (0 or 1) value,
assuming that u represents the onset timings of discrete events in thebrain. To extract discrete events in the brain, our assumption is considered better (Fig. 5D).
Recently, we proposed a general method to estimate MEG/EEG
waveforms that are common across trials (Common Waveform
Estimation, CWE) (Takeda et al., 2010, 2014). CWE can also estimate
spatiotemporal patterns even if their onsets are unknown. STeP is different from CWE mainly in the assumed data structure. CWE is used fortrial data (channel × time × trial), while STeP is used for continuous
data (channel × time). Because resting-state data are continuous, CWE
cannot be used.
Using clustering methods and hidden Markov models, resting-state
MEG/EEG data have been segmented into stationary states (e.g. microstates) (Baker et al., 2014; Britz et al., 2010; Van de Ville et al., 2010;
Woolrich et al., 2013). The differences between STeP and these approaches are particularly evident when propagations occur. If brain activities propagate across different regions (e.g. from V1 to MT), the
MEG/EEG spatial pattern consecutively changes with time. In this case,
STeP represents it with a spatiotemporal pattern, while these approaches represent it by a transition across states. Therefore, STeP is
considered suitable for capturing propagations of brain activities.
Spatiotemporal patterns in resting-state fMRI data
Using STeP, we estimated the two spatiotemporal patterns from theresting-state fMRI data (Fig. 12A). The spatiotemporal patterns showthat the large activities propagate from the auditory cortex to veins
via the visual and cerebellum cortices (Fig. 12A). Recently, Tong et al.
(2012, 2015) showed that fMRI data are strongly contaminated withsLFOs, which are observed at peripheral sites (e.g. ﬁngertip). They also
showed that sLFOs reach different brain regions with different delays(Tong et al., 2012, 2015) and the delays are partially responsible forthe emergence of RSNs (Tong et al., 2015). It is believed that the delays
reﬂect the cerebral blood ﬂow formed by vascular anatomy. The propagation pattern in the estimated spatiotemporal patterns (Fig. 12A)
seems to be consistent with the delay map of sLFOs (Fig. 3A in Tonget al., 2015). Therefore, we suggest that the spatiotemporal patterns
represent how fMRI activities propagate along the cerebral blood ﬂow,
while their origins, whether neuronal or non-neuronal, are not clear.
Applicability
Thinking of an issue, remembering the past, and imagining the future: These activities occur spontaneously and are important brain functions for human beings. However, their details are still veiled because
they occur silently in our brain without exhibiting overt signs. STeP

Y. Takeda et al. / NeuroImage 133 (2016) 251–265

can detect such events from various types of resting-state data, so it
could help us to unveil their mechanisms.
Moreover, applying STeP to continuous data during “tasks” may also
bring interesting results. Usually, continuous MEG/EEG data during
tasks are segmented into trials, and the data during inter-trial intervals
are discarded. However, it is possible that important brain activities also
occur during inter-trial intervals. If so, applying STeP could bring us important ﬁndings. Indeed, we applied STeP to the continuous MEG data
during the visual stimulation task and found that a spatiotemporal pattern resembling the VEF appears not only after the stimulus onsets but
also during the inter-trial intervals (Fig. 9E).
A MATLAB implementation of the proposed method is available
from http://www.cns.atr.jp/takeda/STeP.html.
Acknowledgments
The authors thank Dr. Aapo Hyvärinen for his helpful comments.
This research was supported by contracts with the National Instituteof Information and Communications Technology entitled, “Multimodal
integration for brain imaging measurements” (13601) and “Development of network dynamics modeling methods for human brain data
simulation systems” (173), and JSPS KAKENHI Grant Number 23700499.
Appendix A. Supplementary data
Supplementary data to this article can be found online at http://dx.
doi.org/10.1016/j.neuroimage.2016.03.014.
References
Akaike, H., 1974. A new look at the statistical model identiﬁcation. IEEE Trans. Autom.
Control AC-19, 716–723.
Anemüller, J., Sejnowski, T.J., Makeig, S., 2003. Complex independent component analysis
of frequency-domain electroencephalographic data. Neural Netw. 16, 1311–1323.
Baker, A.P., Brookes, M.J., Rezek, I.A., Smith, S.M., Behrens, T., Probert Smith, P.J., Woolrich,M., 2014. Fast transient networks in spontaneous human brain activity. Elife 3, e01867.
Beckmann, C.F., DeLuca, M., Devlin, J.T., Smith, S.M., 2005. Investigations into resting-stateconnectivity using independent component analysis. Philos. Trans. R. Soc. Lond. Ser. B
Biol. Sci. 360, 1001–1013.
Biswal, B., Yetkin, F.Z., Haughton, V.M., Hyde, J.S., 1995. Functional connectivity in the
motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34,
537–541.
Biswal, B.B., Mennes, M., Zuo, X.N., Gohel, S., Kelly, C., Smith, S.M., Beckmann, C.F.,
Adelstein, J.S., Buckner, R.L., Colcombe, S., Dogonowski, A.M., Ernst, M., Fair, D.,
Hampson, M., Hoptman, M.J., Hyde, J.S., Kiviniemi, V.J., Kotter, R., Li, S.J., Lin, C.P.,
Lowe, M.J., Mackay, C., Madden, D.J., Madsen, K.H., Margulies, D.S., Mayberg, H.S.,
McMahon, K., Monk, C.S., Mostofsky, S.H., Nagel, B.J., Pekar, J.J., Peltier, S.J., Petersen,
S.E., Riedl, V., Rombouts, S.A., Rypma, B., Schlaggar, B.L., Schmidt, S., Seidler, R.D.,
Siegle, G.J., Sorg, C., Teng, G.J., Veijola, J., Villringer, A., Walter, M., Wang, L., Weng,
X.C., Whitﬁeld-Gabrieli, S., Williamson, P., Windischberger, C., Zang, Y.F., Zhang,
H.Y., Castellanos, F.X., Milham, M.P., 2010. Toward discovery science of humanbrain function. Proc. Natl. Acad. Sci. U. S. A. 107, 4734–4739.
Britz, J., Van De Ville, D., Michel, C.M., 2010. BOLD correlates of EEG topography reveal
rapid resting-state network dynamics. NeuroImage 52, 1162–1170.
Brookes, M.J., Woolrich, M., Luckhoo, H., Price, D., Hale, J.R., Stephenson, M.C., Barnes, G.R.,
Smith, S.M., Morris, P.G., 2011. Investigating the electrophysiological basis of resting265

state networks using magnetoencephalography. Proc. Natl. Acad. Sci. U. S. A. 108,
16783–16788.
de Cheveigné, A., Simon, J.Z., 2007. Denoising based on time-shift PCA. J. Neurosci.
Methods 165, 297–305.
de Pasquale, F., Della Penna, S., Snyder, A.Z., Lewis, C., Mantini, D., Marzetti, L., Belardinelli,P., Ciancetta, L., Pizzella, V., Romani, G.L., Corbetta, M., 2010. Temporal dynamics of
spontaneous MEG activity in brain networks. Proc. Natl. Acad. Sci. U. S. A. 107,
6040–6045.
de Pasquale, F., Della Penna, S., Snyder, A.Z., Marzetti, L., Pizzella, V., Romani, G.L., Corbetta,
M., 2012. A cortical core for dynamic integration of functional networks in the resting
human brain. Neuron 74, 753–764.
Delorme, A., Makeig, S., 2004. EEGLAB: an open source toolbox for analysis of single-trial
EEG dynamics including independent component analysis. J. Neurosci. Methods 134,
9–21.
Deuker, L., Olligs, J., Fell, J., Kranz, T.A., Mormann, F., Montag, C., Reuter, M., Elger, C.E.,
Axmacher, N., 2013. Memory consolidation by replay of stimulus-speciﬁc neuralactivity. J. Neurosci. 33, 19373–19383.
Dragoi, G., Tonegawa, S., 2011. Preplay of future place cell sequences by hippocampal cellular assemblies. Nature 469, 397–401.
Dragoi, G., Tonegawa, S., 2013. Distinct preplay of multiple novel spatial experiences in
the rat. Proc. Natl. Acad. Sci. U. S. A. 110, 9100–9105.
Foster, D.J., Wilson, M.A., 2006. Reverse replay of behavioural sequences in hippocampal
place cells during the awake state. Nature 440, 680–683.
Fox, M.D., Raichle, M.E., 2007. Spontaneous ﬂuctuations in brain activity observed withfunctional magnetic resonance imaging. Nat. Rev. Neurosci. 8, 700–711.
Fox, M.D., Snyder, A.Z., Vincent, J.L., Corbetta, M., Van Essen, D.C., Raichle, M.E., 2005. The
human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proc. Natl. Acad. Sci. U. S. A. 102, 9673–9678.
Han, F., Caporale, N., Dan, Y., 2008. Reverberation of recent visual experience in spontaneous cortical waves. Neuron 60, 321–327.
Hopﬁeld, J.J., 1995. Pattern recognition computation using action potential timing forstimulus representation. Nature 376, 33–36.
Ikegaya, Y., Aaron, G., Cossart, R., Aronov, D., Lampl, I., Ferster, D., Yuste, R., 2004. Synﬁrechains and cortical songs: temporal modules of cortical activity. Science 304, 559–564.
Izhikevich, E.M., 2006. Polychronization: computation with spikes. Neural Comput. 18,
245–282.
Ji, D., Wilson, M.A., 2007. Coordinated memory replay in the visual cortex and hippocampus during sleep. Nat. Neurosci. 10, 100–107.
Jung, T.P., Makeig, S., Westerﬁeld, M., Townsend, J., Courchesne, E., Sejnowski, T.J., 2001.
Analysis and visualization of single-trial event-related potentials. Hum. Brain Mapp.
14, 166–185.
MacQueen, J., 1967. Some methods for classiﬁcation and analysis of multivariate observations. In: LeCam, L.M., Neyman, J. (Eds.), Proceedings of the Fifth Berkeley Symposiumon Mathematical Statistics and Probability 1. University of California Press,
pp. 281–297.
Majeed, W., Magnuson, M., Hasenkamp, W., Schwarb, H., Schumacher, E.H., Barsalou, L.,
Keilholz, S.D., 2011. Spatiotemporal dynamics of low frequency BOLD ﬂuctuations
in rats and humans. NeuroImage 54, 1140–1150.
Mantini, D., Perrucci, M.G., Del Gratta, C., Romani, G.L., Corbetta, M., 2007. Electrophysiological signatures of resting state networks in the human brain. Proc. Natl. Acad. Sci.
U. S. A. 104, 13170–13175.
Olshausen, B.A., 2003. Learning sparse, overcomplete representations of time-varying
natural images. IEEE International Conference on Image Processing.
Raichle, M.E., MacLeod, A.M., Snyder, A.Z., Powers, W.J., Gusnard, D.A., Shulman, G.L.,
2001. A default mode of brain function. Proc. Natl. Acad. Sci. U. S. A. 98, 676–682.
Schwarz, G., 1978. Estimating the dimension of a model. Ann. Stat. 6, 461–464.
Smith, S.M., Fox, P.T., Miller, K.L., Glahn, D.C., Fox, P.M., Mackay, C.E., Filippini, N., Watkins,
K.E., Toro, R., Laird, A.R., Beckmann, C.F., 2009. Correspondence of the brain's functionalarchitecture during activation and rest. Proc. Natl. Acad. Sci. U. S. A. 106, 13040–13045.
Staresina, B.P., Alink, A., Kriegeskorte, N., Henson, R.N., 2013. Awake reactivation predicts
memory in humans. Proc. Natl. Acad. Sci. U. S. A. 110, 21159–21164.
Takeda, Y., Sato, M., Yamanaka, K., Nozaki, D., Yamamoto, Y., 2010. A generalized method
to estimate waveforms common across trials from EEGs. NeuroImage 51, 629–641.
Takeda, Y., Yamanaka, K., Yamagishi, N., Sato, M., 2014. Revealing time-unlocked brain activity from MEG measurements by common waveform estimation. PLoS One 9, e98014.
Tambini, A., Ketz, N., Davachi, L., 2010. Enhanced brain correlations during rest are related
to memory for recent experiences. Neuron 65, 280–290.
Tong, Y., Hocke, L.M., Licata, S.C., Frederick, B., 2012. Low-frequency oscillations measured
in the periphery with near-infrared spectroscopy are strongly correlated with blood
oxygen level-dependent functional magnetic resonance imaging signals. J. Biomed.
Opt. 17, 106004.
Tong, Y., Hocke, L.M., Fan, X., Janes, A.C., Frederick, B., 2015. Can apparent resting stateconnectivity arise from systemic ﬂuctuations? Front. Hum. Neurosci. 9, 285.
Van de Ville, D., Britz, J., Michel, C.M., 2010. EEG microstate sequences in healthy humans
at rest reveal scale-free dynamics. Proc. Natl. Acad. Sci. U. S. A. 107, 18179–18184.
Wilson, M.A., McNaughton, B.L., 1994. Reactivation of hippocampal ensemble memories
during sleep. Science 265, 676–679.
Woolrich, M.W., Baker, A., Luckhoo, H., Mohseni, H., Barnes, G., Brookes, M., Rezek, I.,
2013. Dynamic state allocation for MEG source reconstruction. NeuroImage 77,
77–92.