# State-transition dynamics of resting-state functional magnetic resonance
imaging data: Model comparison and test-to-retest analysis

Saiful Islam1, Pitambar Khanra2, Johan Nakuci4, Sarah F. Muldoon1,2,3,
Takamitsu Watanabe5, Naoki Masuda1,2∗ 1Institute for Artificial Intelligence
and Data Science, University at Buffalo, State University of New York at
Buffalo, USA
2Department of Mathematics, University at Buffalo, State University of New
York at Buffalo, USA
3Neuroscience Program, University at Buffalo, State University of New York at
Buffalo, USA
4School of Psychology, Georgia Institute of Technology, USA
5International Research Centre for Neurointelligence, The University of Tokyo
Institutes for Advanced Study, Tokyo, Japan
∗Corresponding author<EMAIL_ADDRESS>
###### Abstract

Electroencephalogram (EEG) microstate analysis entails finding dynamics of
quasi-stable and generally recurrent discrete states in multichannel EEG time
series data and relating properties of the estimated state-transition dynamics
to observables such as cognition and behavior. While microstate analysis has
been widely employed to analyze EEG data, its use remains less prevalent in
functional magnetic resonance imaging (fMRI) data, largely due to the slower
timescale of such data. In the present study, we extend various data
clustering methods used in EEG microstate analysis to resting-state fMRI data
from healthy humans to extract their state-transition dynamics. We show that
the quality of clustering is on par with that for various microstate analyses
of EEG data. We then develop a method for examining test-retest reliability of
the discrete-state transition dynamics between fMRI sessions and show that the
within-participant test-retest reliability is higher than between-participant
test-retest reliability for different indices of state-transition dynamics,
different networks, and different data sets. This result suggests that state-
transition dynamics analysis of fMRI data could discriminate between different
individuals and is a promising tool for performing fingerprinting analysis of
individuals.

###### keywords:

fMRI, EEG, MEG, microstates, clustering, dynamics, test-retest reliability,
fingerprinting.

##  1 Introduction

Activity of the human brain is dynamic even at rest, and brain dynamics on
various spatial scales are considered to drive myriad functions of the brain
[1, 2, 3, 4]. Multiple methods to characterize brain dynamics have been
proposed, many of which rely on the detection of brain states and
quantification of how the brain transitions through such states. Microstate
analysis is an early-proposed method for estimating discrete states in
electroencephalogram (EEG) data [5, 6, 7]. EEG microstate analysis usually
entails clustering of multi-electrode EEG signals, with each data point to be
clustered corresponding to a time point of the measurement. Each cluster, or
microstate, is a representation of a global functional state of the brain.
Microstates obtained from resting-state EEG data tend to last about 100 ms and
are reproducible [8, 9, 6, 10, 11]. Microstate analysis has been extended for
magnetoencephalography (MEG) data, with the microstates being estimated by
conventional clustering methods [12, 13] or the hidden-Markov model (HMM) [14,
12] among other methods.

Microstate analysis in its original sense (i.e., detecting and utilizing
microstates lasting about 100 ms) does not directly apply to functional
magnetic resonance imaging (fMRI) data because the temporal resolution of fMRI
is limited, preventing one from detecting dynamics on the timescale of 100 ms.
One direction to resolve this limitation is to use EEG microstate analysis
results to inform states in fMRI data [15, 16, 17, 18]. An alternative
approach is to estimate and use state-transition dynamics of spatial fMRI
signals, as microstate analysis does for EEG (and MEG) data, regardless of
different time resolutions between fMRI and EEG/MEG. Such state-transition
dynamics for fMRI data have been estimated by data clustering algorithms as in
the case of the EEG/MEG microstate analysis [19, 20, 21, 22], the HMM or its
variants [23, 24, 25, 26, 21, 27, 22], and energy landscape analysis [28, 29,
30, 31]. Each discrete state in fMRI data corresponds to a vector of activity
patterns at specified regions of interests (ROIs) [32, 33, 34, 22, 35], or a
functional network among ROIs [36, 37, 38, 20, 19, 39, 23, 25, 21].

In general, successful single individual inferences from neuroimaging data
would suggest their potential applications for both scientific investigations
and clinical practice. Research has shown that functional networks from fMRI
data can be used as a reliable fingerprint of human individuals through test-
retest analyses [40, 41, 42, 43, 44, 45]. Test-retest reliability has also
been assessed for dynamic functional networks estimated from fMRI data [46,
47, 48], whereas test-retest reliability for dynamic functional networks has
been reported to be lower than that for static functional networks [46, 47].
With this study, we are interested in test-retest reliability of state-
transition dynamics in fMRI data, which has been underexplored.

In the present study, we assess the potential effectiveness of dynamics of
discrete states estimated from fMRI data at fingerprinting individuals. Here,
we use fMRI data as multivariate time series, each dimension of which
represents a single ROI, akin to microstate analysis for EEG and MEG data.
This approach contrasts with the aforementioned prior studies on test-retest
reliability of dynamic functional networks. Our analysis involves examination
of what methodological choices (e.g., the clustering method applied to the
fMRI data to define discrete states, the number of clusters identified, and
the indices used to characterize the estimated state transition dynamics)
yield a higher test-retest reliability of the state-transition dynamics; such
an assessment has previously been carried out for EEG microstate analysis [9].
Based on a permutation test to quantify test-retest reliability, we show that,
in general, transitory dynamics of discrete states estimated for fMRI data
yield higher within-participant than between-participant test-retest
reliability across clustering methods, the number of clusters, observables of
the state-transition dynamics, two sets of ROIs, and two data sets. Code for
computing dynamics of discrete states and their test-retest reliability used
in the present paper is available on Github [49].

##  2 Methods

###  2.1 Midnight Scan Club data

We use the resting-state fMRI data provided by the Midnight Scan Club (MSC)
project [41]. The MSC’s resting-state fMRI data consist of recording from ten
healthy human adults over ten consecutive nights. A single recording session
of the resting-state fMRI experiment lasted for 30 mins, resulting in 818
volumes. The imaging was performed on a Siemens TRIO 3T MRI scanner. All
functional imaging was performed using an echo planar imaging (EPI) sequence
(TR $=2.2$ s, TE $=27$ ms, flip angle $=90^{\circ}$, voxel size $=4$ mm
$\times$ 4 mm $\times$ 4 mm, 36 slices).

It was originally reported that the eighth participant (i.e., MSC08) fell
asleep, showed frequent and prolonged eye closures, and had systematically
large head motion, yielding considerably less reliable data than those
obtained from the other participants [41]. In our previous work, we also
noticed that the quality of the data analysis fluctuated considerably more
across the different sessions for the tenth participant (i.e., MSC10) than for
the other participants except MSC08 [50]. Therefore, we excluded MSC08 and
MSC10 from the following analysis.

We used SPM12 (http://www.fil.ion.ucl.ac.uk/spm) to preprocess the resting-
state functional images. Specifically, we first conducted realignment,
unwraping, slice-timing correction, and normalization to a standard template
(ICBM 152). Then, we performed regression analyses to remove the effects of
head motion, white matter signals, and cerebrospinal fluid signals. Lastly, we
conducted band-pass temporal filtering (0.01–0.1 Hz).

We used a DMN composed of 12 ROIs [51]. To optionally reduce the dimension of
the DMN, we averaged over each pair of the symmetrically located right- and
left-hemisphere ROIs into one observable. The symmetrized DMN has eight ROIs
because four ROIs (i.e., amPFC, vmPFC, pCC, and retro splen) in the original
coordinate system are approximately on the midline and therefore have not
undergone the averaging over the right- and left-hemisphere ROIs [31].

In addition to the DMN, we also analyzed the so-called whole-brain network. We
determined the regions of interest (ROIs) of the whole-brain network by
employing the 264 spherical ROIs whose coordinates were identified in a
previous study [52]. We then removed 50 ROIs labelled ‘uncertain’ or
‘subcortical’, resulting in 214 ROIs. The 214 ROIs were labeled either of the
following nine functionally different brain systems: auditory network, dorsal
attention network (DAN), ventral attention network (VAN), cingulo-opercular
network (CON), default mode network (DMN), fronto-parietal network (FPN),
salience network (SAN), somatosensory and motor network (SMN), or visual
network. We merged the DAN, VAN, and CON into an attention network (ATN) to
reduce the number of observables from nine to seven, as we did in our previous
studies [53, 50]. This is because the DAN, VAN, and CON have been suggested to
be responsible for similar attention-related cognitive activity [52]. We then
calculated the average fMRI signal for each of the seven systems by first
averaging the signal over the volumes in the sphere of radius 4 mm centered
around the provided coordinate of each ROI [52], and then averaging the signal
over all ROIs belonging to the system (e.g., 13 ROIs in the auditory network).
We call the thus obtained seven-dimensional system the whole-brain network. It
should be noted that the DMN constitutes one ROI in the whole-brain network,
whereas the DMN described above as a system of ROIs is composed of either 8 or
12 ROIs depending on whether or not we average over symmetrically located
ROIs.

###  2.2 Human Connectome Project data

We also analyzed the fMRI data recorded from healthy human adults and shared
as the S1200 data in the Human Connectome Project (HCP) [54]. In the S1200
data, 1200 adults between $22$–$35$ years old underwent four sessions of
15-min EPI sequence with a 3T Siemens Connectome-Skyra (TR $=0.72$ s, TE
$=33.1$ ms, $72$ slices, $2.0$ mm isotropic, field of view (FOV) $=208\times
180$ mm) and a T1-weighted sequence (TR $=2.4$ s, TE $=2.14$ ms, $0.7$ mm
isotropic, FOV $=224\times 224$ mm). Here, we only analyzed the 100 unrelated
participant subset released by the HCP. All these 100 participants completed
both diffusion weighted MRI and two resting-state fMRI scans.

Each participant underwent two sessions of resting-state fMRI recording, and
each session consisted of both Left-Right (LR) and Right-Left (RL) phases. In
the following text, we refer to phases as sessions. Therefore, each
participant’s data consist of four sessions. We used data from participants
with at least 1150 volumes in each of the four sessions after we had removed
volumes with motion artifacts, resulting in a final analysis of 87
participants. For the 87 participants, we removed the volumes with motion
artifacts and then used the last 1150 volumes in each session, with the aim of
removing possible transient effects.

We employed independent component analysis (ICA) to remove nuisance and motion
signals [55]. Then, any volumes with frame displacement greater than 0.2 mm
[56] were excised [57]. This is because the ICA-FIX pipeline has been found
not to fully remove motion-related artifacts [58, 59]. Next, we standardized
each voxel by subtracting the temporal mean, and then global signal regression
(see section 2.3) was carried out.

We averaged the fMRI signal over all the voxels within each ROI of the AAL
atlas [60] in each volume. We remark that the AAL atlas is composed of 116
ROIs. In order to map these ROIs to representative brain systems, we first
mapped each of the cortical ROIs to the parcellation scheme from the
Schaefer-100 atlas [61]. We based the assignment of the ROI to the brain
system that minimized the Euclidian distance from the centroid of an ROI in
the AAL to the corresponding centroid of an ROI in the Schaefer atlas. After
we assigned each ROI to a system, we removed 42 ROIs labeled ‘subcortical’ or
‘cerebellar’, which yielded 74 ROIs. These 74 ROIs were then assigned to one
of the $N=7$ functionally different brain networks: control network, DMN, DAN,
limbic network, salience/ventral attention network, somatomotor network, and
visual network. We call this seven-dimensional system the whole-brain network
for the HCP data. Similarly to the case of the whole-brain network for the MSC
data, we first averaged the fMRI signal over the voxels within each ROI and
then further averaged the signal over the ROIs belonging to the same system
(e.g., 59 ROIs belonging to the DMN).

###  2.3 Global signal removal

We denote the fMRI time series for a session by
$\\{\bm{x}_{1},\ldots,\bm{x}_{T}\\}$, where $T$ is the number of volumes
(i.e., time points), $\bm{x}_{t}=(x_{t,1},\ldots,x_{t,\tilde{N}})$ is the fMRI
signal at time $t$, and $\tilde{N}$ is the number of ROIs with which we
compute the global signal. Note that $\tilde{N}$ may be larger than $N$, which
occurs when we define a global signal widely from the brain including ROIs
that we do not use for estimating the discrete states. The global signal is
the average of the signal over all the $\tilde{N}$ ROIs at each time, i.e.,

$\overline{x}_{t}=\frac{\sum_{i=1}^{\tilde{N}}x_{t,i}}{\tilde{N}}.$ (1)

We remove the global signal [62] by subtracting $\overline{x}_{t}$ from each
$x_{t,i}$ (with $i\in\\{1,\ldots,\tilde{N}\\}$) and dividing the result by the
standard deviation, i.e.,

$\sigma_{t}=\sqrt{\frac{\sum_{i=1}^{\tilde{N}}(x_{t,i}-\overline{x}_{t})^{2}}{\tilde{N}}}.$
(2)

We carry out this procedure for each $t$.

The global signal in resting-state fMRI data is considered to primarily
consist of physiological noise stemming from different factors such as
respiration, scanner-related artifacts, and motion-related artifacts. By
removing the global signal, several quality-control metrics are improved, the
anatomical specificity of functional-connectivity patterns is enhanced, and
there is a potential increase in behavioral variance [63, 64].

For the DMN obtained from the MSC data, we first removed the global signal
calculated over the $\tilde{N}=30$ ROIs in the coordinate system provided by
[51], which included the $N=12$ ROIs in the DMN. Then, we compared three
treatments of global signal removal for the DMN as follows. In the first and
second treatments, we then removed the global signal calculated from the
$\tilde{N}=12$ DMN ROIs from each of the $12$ ROIs in the DMN. Next, we
averaged the obtained time series over each symmetric pair of DMN ROIs
corresponding to the two hemispheres. If the ROI is roughly on the midline,
there is no such symmetric pair of ROIs, in which case we only removed the
global signal. After aggregating the symmetric pairs of ROIs in this manner,
there are $N=8$ ROIs in the DMN. This concludes the first treatment. In the
case of the second treatment, we additionally removed the global signal
calculated over the $\tilde{N}=N=8$ ROIs. In the third treatment, after the
removal of the global signal calculated over $30$ ROIs, which is common for
all the three treatments, we further removed global signal calculated from the
$\tilde{N}=12$ DMN ROIs from each of the 12 ROIs. We do not further process
the data. Therefore, with the third treatment, the final DMN consists of
$N=12$ ROIs.

For the whole-brain network obtained from the MSC data, we first removed the
global signal computed from the $\tilde{N}=264$ ROIs. Then, we extracted $N=7$
-dimensional time series as described in section 2.1. Finally, we further
removed the global signal computed from the $\tilde{N}=N=7$ ROIs in the whole-
brain network. The global signal removal for the whole-brain network obtained
from the HCP data is the same except that we computed the first global signal
from the $\tilde{N}=116$ ROIs of the AAL atlas (see section 2.2).

###  2.4 Estimation of discrete states

There are various methods for estimating microstates in the EEG and MEG data
[5, 7, 6, 13]. We tailor seven popular methods for finding microstates in EEG
and MEG data to the case of fMRI data to estimate their discrete states.
Because the discrete states that we find for fMRI data are not equivalent to
EEG/MEG microstates, we refer to the former as states, discrete states, or
clusters in the following text. We describe each method in the following
subsections. See Table 1 for main notations.

Symbols |  Description
---|---
$N_{\text{p}}$ |  Number of participants
$N_{\text{s}}$ |  Number of sessions for each participant
$N$ |  Number of ROIs
$T$ |  Number of volumes (i.e., time points) in each session
$\bm{x}_{t}\in\mathbb{R}^{N}$ |  fMRI signal at time $t$
$\overline{x}_{t}$ |  Average of $x_{t,i}$ over ROIs
$\sigma_{t}$ |  Standard deviation of $x_{t,i}$ over ROIs
$K$ |  Number of discrete states
$\bm{c}_{\ell}\in\mathbb{R}^{N}$ |  Centroid of the $\ell$th cluster, where $\ell\in\\{1,2,\ldots,K\\}$
$L_{t}$ |  Cluster label for $\bm{x}_{t};L_{t}\in\\{1,2,\ldots,K\\}$
Table 1: Main notations used in this paper.

####  2.4.1 K-means clustering

The K-means clustering is a simple and popular clustering method to partition
the data points into $K$ mutually exclusive clusters. Various EEG and MEG
microstate analysis [6, 65, 13] and the studies on temporal variability of
functional connectivity states of fMRI data [39, 66, 67] used the K-means
clustering. It starts with a predefined number of clusters, $K$. We initialize
the centroids of the clusters by the k-means++ algorithm [68]. The k-means++
algorithm consists of the following steps. In step (i), we select one centroid
uniformly at random from all the data points. In step (ii), for each data
point $\bm{x}_{t}$ that is not yet selected as a centroid, we calculate the
distance to the nearest centroid. In step (iii), we sample one $\bm{x}_{t}$
that is not yet selected as a centroid yet with the probability proportional
to the square of the distance between $\bm{x}_{t}$ to the nearest centroid. In
step (iv), we add the $\bm{x}_{t}$ sampled in the third step as a new
centroid. We repeat steps (ii), (iii), and (iv) until we obtain $K$ centroids.
This initialization method accelerates the convergence of the algorithm.

Then, we refine the $K$ centroids, denoted by $\bm{c}_{1}$, $\ldots$,
$\bm{c}_{K}$, as follows. The first step is to assign each data point
$\bm{x}_{t}$ to the nearest centroid, i.e., the centroid realizing

$L_{t}=\text{arg}\min_{\ell\in\\{1,\ldots,K\\}}\left\|\bm{x}_{t}-\bm{c}_{\ell}\right\|,$
(3)

where $\left\|\cdot\right\|$ denotes the Euclidean norm. The second step is to
update the centroid of each cluster $\ell$ by the average of the data points
belonging to the cluster as follows:

$\bm{c}_{\ell}=\frac{\sum_{t=1}^{T}\bm{x}_{t}\delta_{L_{t},\ell}}{\sum_{t=1}^{T}\delta_{L_{t},\ell}},$
(4)

where $\delta_{L_{t},\ell}$ is the Kronecker delta; $\delta_{L_{t},\ell}=1$ if
$L_{t}=\ell$, and $\delta_{L_{t},\ell}=0$ otherwise. The Kronecker delta in
the equation allows us to take the summation only over the data points
belonging to the $\ell$th cluster. We repeat the first and second steps until
the change in the residual sum of squares (RSS), defined by

$\text{RSS}=\sum_{t=1}^{T}\sum_{\ell=1}^{K}\delta_{L_{t},\ell}\left\|\bm{x}_{t}-\bm{c}_{\ell}\right\|^{2},$
(5)

between the two subsequent steps falls below $10^{-5}$ for the first time. We
use the implementation of $k$-means in scikit-learn [69].

####  2.4.2 K-medoids clustering

The K-medoids clustering algorithm [70] is a variant of the K-means
clustering. The K-medoids clustering uses the original data points as the
centroids of the clusters, referred to as medoids. In contrast, the K-means
clustering uses the average of the points in the cluster as the centroid of
the cluster. The K-medoids clustering begins with a set of $K$ data points as
medoids, which we select using the k-medoids++ method. In fact, k-medoids++ is
the same as k-means++. In the next step, we assign each $\bm{x}_{t}$ to the
$\ell$th cluster whose medoid is closest to $\bm{x}_{t}$ in terms of the
Euclidean distance. Then, we update the medoid of each cluster to $\bm{x}_{t}$
that belongs to the cluster and minimizes the sum of the Euclidean distance to
the other data points in the same cluster. We repeat the last two steps until
the dissimilarity score (i.e., the sum of the Euclidean distance from the
medoid to the other data points in the cluster) stops changing for each
cluster. We use the $k$-medoids implemented in scikit-learn [69].

####  2.4.3 Agglomerative hierarchical clustering

Agglomerative hierarchical clustering, which we simply call the agglomerative
clustering (AC), is a bottom-up clustering method. The AC method initially
regards each data point as a single-node cluster. Then, one merges a pair of
clusters one after another based on a linkage criterion. Among various linkage
criteria, we use the Ward’s method implemented in scikit-learn [69]. In each
step of merging two clusters, the Ward’s method minimizes the within-cluster
variance, i.e., the squared Euclidean distance between $\bm{x}_{t}$ and the
centroid of the new cluster to which $\bm{x}_{t}$ belongs, which is summed
over $t\in\\{1,\ldots,T\\}$. We stop the cluster merging procedure once the
number of clusters is equal to $K$.

####  2.4.4 Atomize and agglomerate hierarchical clustering

The Atomize and Agglomerate Hierarchical Clustering (AAHC) is another bottom-
up hierarchical clustering algorithm [71, 72, 9]. A main difference between
AAHC and traditional bottom-up hierarchical clustering methods is that AAHC
atomizes the worst cluster. In other words, AAHC disintegrates the worst
cluster and assigns each member of this cluster to a different cluster instead
of merging the entire worst cluster with the most similar cluster.

AAHC uses the global explained variance (GEV) as a measure of the quality of
the cluster [72, 9, 6, 73]. The GEV for the $\ell$th cluster is defined by

$\text{GEV}_{\ell}=\frac{\sum_{t=1}^{T}\delta_{L_{t},\ell}\,\text{corr}(\bm{x}_{t},\bm{c}_{\ell})^{2}\,\sigma_{t}^{2}}{\sum_{t=1}^{T}\sigma_{t}^{2}},$
(6)

where $\text{corr}(\bm{x}_{t},\bm{c}_{\ell})$ is the cosine similarity between
$\bm{x}_{t}$ and $\bm{c}_{t}$ given by

$\text{corr}(\bm{x}_{t},\bm{c}_{\ell})=\frac{\langle\bm{x}_{t},\bm{c}_{\ell}\rangle}{\left\|\bm{x}_{t}\right\|\left\|\bm{c}_{\ell}\right\|}.$
(7)

In Eq. (7), $\langle\bm{x}_{t},\bm{c}_{\ell}\rangle$ is the inner product of
$\bm{x}_{t}$ and $\bm{c}_{\ell}$. Variable $\sigma_{t}$ represents the
standard deviation of the data point $\bm{x}_{t}$ across the ROIs and is given
by Eq. (2). Quantity $\sigma_{t}$ is known as global field power (GFP) in the
literature of microstate analysis for EEG and MEG data [74, 72, 5, 13]. For
the second and third treatments of the global signal removal, it holds true
that $\sigma_{t}=1$ for any $t$ because of the global signal removal carried
out in the last step of the treatment.

In the AAHC, we define the worst cluster as the one with the smallest
$\text{GEV}_{\ell}$, $\ell\in\\{1,2,\cdots,K\\}$ and atomize it. Then, we
assign each data point $\bm{x}_{t}$ of the atomized cluster to the $\ell$th
cluster that maximizes Eq. (7) [72, 9]. As in the AC, the AAHC initially
regards each $\bm{x}_{t}$ as a single-node cluster. We repeat finding the
worst cluster, atomizing it, and assigning each $\bm{x}_{t}$ in the atomized
cluster to a different cluster until the number of clusters reaches $K$.

####  2.4.5 Topographic atomize and agglomerate hierarchical clustering

The Topographic Atomize and Agglomerate Hierarchical Clustering (TAAHC) is a
modification of AAHC [9, 73]. The difference between AAHC and TAAHC is that
TAAHC defines the worst cluster to be the $\ell$th cluster that is the
smallest in terms of the sum of the correlation of the data points in the
cluster with its centroid $\bm{c}_{\ell}$ [9, 73]. In other words, the worst
cluster $\ell$ is the minimizer of

$\displaystyle\text{CRS}(\ell)=\sum_{t=1}^{T}\delta_{L_{t},\ell}\,\text{corr}(\bm{x}_{t},\bm{c}_{\ell})=\sum_{t=1}^{T}\frac{\delta_{L_{t},\ell}\,\langle\bm{x}_{t},\bm{c}_{\ell}\rangle}{\left\|\bm{x}_{t}\right\|\left\|\bm{c}_{\ell}\right\|}$
(8)

over $\ell\in\\{1,\ldots,K\\}$. As in the AC and AAHC, the TAAHC first regards
each $\bm{x}_{t}$ as a single-node cluster. Second, we identify the cluster
with the smallest $\text{CRS}(\ell)$. Third, we atomize the selected cluster
and reassign each of its member $\bm{x}_{t}$ to the cluster whose centroid is
the closest to $\bm{x}_{t}$ in terms of
$\text{corr}(\bm{x}_{t},\bm{c}_{\ell})$. We iterate the second and third steps
until we obtain $K$ clusters.

####  2.4.6 Bisecting K-means clustering

The bisecting K-means method combines the K-means clustering method and
divisive hierarchical clustering [75]. Initially, we let all data points form
a single cluster. Then, we apply the K-means clustering with $K=2$ to
partition the data points into two clusters, following the procedure described
in section 2.4.1. Then, we select the cluster that has the larger value of the
dissimilarity defined for the $\ell$th cluster by

$\rm{SSE}_{\ell}=\sum_{t=1}^{T}\delta_{L_{t},\ell}\left\|\bm{x}_{t}-\bm{c}_{\ell}\right\|^{2}.$
(9)

Then, we run the K-means clustering on the selected cluster to split it into
two clusters. We repeat selecting the cluster with the largest
$\rm{SSE}_{\ell}$ and bisecting it until we obtain $K$ clusters. We use the
implementation of the bisecting K-means in scikit-learn [69].

####  2.4.7 Gaussian mixture model

The Gaussian mixture models (GMM) represents each cluster as a multivariate
Gaussian distribution. We denote by $N(\bm{\mu}_{\ell},\bm{\Sigma}_{\ell})$,
with $\ell\in\\{1,2,\ldots,K\\}$, the multidimensional Gaussian distribution
with mean vector $\bm{\mu}_{\ell}$ and covariance matrix $\bm{\Sigma}_{\ell}$
[76, 22]. The GMM is given by

$p(\bm{x}_{t})=\sum_{\ell=1}^{K}\pi_{\ell}\,\mathcal{N}(\bm{x}_{t}|\bm{\mu}_{\ell},\bm{\Sigma}_{\ell}),$
(10)

where $\pi_{\ell}$ is the mixing weight, i.e., the probability that a data
point originates from the $\ell$th multivariate Gaussian distribution. Note
that $\sum_{\ell=1}^{K}\pi_{\ell}=1$. The likelihood function for the set of
all the data points is given by

$p(\bm{x}_{1},\ldots,\bm{x}_{T})=\prod_{t=1}^{T}\sum_{\ell=1}^{K}\pi_{\ell}\,\mathcal{N}(\bm{x}_{t}|\bm{\mu}_{\ell},\bm{\Sigma}_{\ell}).$
(11)

We infer the parameter values by maximizing the log-likelihood function using
an expectation-maximization (EM) algorithm [77, 76, 22]. We regard
$\bm{\mu}_{\ell}$ as the centroid of the $\ell$th cluster. Because the GMM is
a soft clustering method, we assign each time point $t$ to the $\ell$th
cluster that maximizes
$\hat{\pi}_{\ell}\mathcal{N}(\bm{x}_{t}|\hat{\bm{\mu}}_{\ell},\hat{\bm{\Sigma}}_{\ell})$,
where $\hat{\pi}_{\ell}$, $\hat{\bm{\mu}}_{\ell}$, and
$\hat{\bm{\Sigma}}_{\ell}$ are the obtained maximum likelihood estimator. We
use the GaussianMixture class in scikit-learn, which uses K-means clustering
for initializing the parameters [69].

Among the seven methods that we employ to cluster the fMRI data, the GMM is
the only parametric model. All the other methods are non-parametric clustering
methods.

###  2.5 Evaluation of the clustering methods

The number of microstates estimated for EEG and MEG data depends on studies
[74, 9, 6, 65]. Studies on temporal dynamics of functional connectivity in
fMRI data are also diverse with respect to the number of clusters [39, 66,
67]. Therefore, we examine the number of states, $K$, from 2 to 10 for each
clustering algorithm. To compare the quality of the different clustering
methods, we use the GEV given by Eq. (6). The GEV captures the amount of the
data variance explained by the microstates’ centroids, also called the global
map, cluster map, microstate map, and template map [72, 9, 7, 13]. We
calculate the total GEV as the sum of the GEV over all the states, i.e.,

$\text{GEV}_{\text{total}}=\sum_{\ell=1}^{K}\text{GEV}_{\ell}$ (12)

and average it over all the sessions and participants. A large value of the
$\text{GEV}_{\text{total}}$ suggests that the obtained clustering is of high
quality.

We also measure the quality of the clustering methods using the within-cluster
sum of squares (WCSS) [78], also known as the distortion measure [76]. The
WCSS is defined by

$\text{WCSS}=\sum_{t=1}^{T}\sum_{\ell=1}^{K}\delta_{L_{t},\ell}\left\|\bm{x}_{t}-\bm{c}_{\ell}\right\|^{2}.$
(13)

A small WCSS value indicates that the data points are tightly clustered and
therefore the clustering is of high quality.

###  2.6 Comparison of state-transition dynamics between different sessions

####  2.6.1 Observables for the state-transition dynamics

To test reproducibility of the fMRI state-transition dynamics across
participants and sessions, we measure the following five observables for each
session. These observables are often used in the analysis of microstate
dynamics for EEG and MEG data [14, 9, 65, 13] and activity patterns for fMRI
data [32, 39, 66].

First, we use the centroid of each of the $K$ states as an observable. The
centroid $\bm{c}_{\ell}$ of the $\ell$th state represents the set of data
points which are assigned to the $\ell$th state. We remind that the centroid
is an $N$-dimensional vector.

Second, the coverage time of the $\ell$th state is the number of times
$t\in\\{1,\ldots,T\\}$ in which the $\ell$th state appears. We normalize the
coverage time of each state by dividing it by the total observation time, $T$.

Third, we measure the frequency of appearance of each state. If the $\ell$th
state starts and then lasts for some time steps before transiting to a
different state, then we say that this is a unique appearance of $\ell$. That
is, we count the consecutive appearances as one unique appearance. The
frequency of appearance of $\ell$ is defined as the number of unique
appearance divided by $T$.

Fourth, the average lifespan of the $\ell$th state is the time spent in a
unique appearance of $\ell$ that is averaged over all unique appearances of
$\ell$. The average lifespan of $\ell$ is equal to the coverage time divided
by the number of unique appearance of $\ell$.

Fifth, we investigate the frequency of transitions from one state to another
as follows. Let $n_{\ell\ell^{\prime}}$ be the number of times with which the
transition from the $\ell$th state to the $\ell^{\prime}$th state occurs in
the given session, where $\ell^{\prime}\neq\ell$. We define the transition
probability from $\ell$ to $\ell^{\prime}$ by
$p_{\ell\ell^{\prime}}=n_{\ell\ell^{\prime}}/\sum_{\ell^{\prime\prime}=1;\ell^{\prime\prime}\neq\ell}^{K}n_{\ell\ell^{\prime\prime}}$
and we set $p_{\ell\ell}=0$. The $K\times K$ transition probability matrix is
given by $P=(p_{\ell\ell^{\prime}})$ with
$\ell,\ell^{\prime}\in\\{1,\ldots,K\\}$.

####  2.6.2 Discrepancy measures for comparing the state-transition dynamics
between two sessions

For examining the reproducibility of state-transition dynamics between
sessions of the same participant and between different participants, we need
to compare observables between pairs of sessions. To this end, we first need
to find the best matching of the states between the two sessions. For
$K\in\\{2,\ldots,8\\}$, we assess all the $K!$ pairwise matchings of the
states between the two sessions. For each matching, we calculate the
correlation between centroids $\bm{c}_{\ell}$ and $\bm{c}_{\ell^{\prime}}$ of
the matched states, i.e., $\ell$th state in the first session and the
$\ell^{\prime}$th state in the second session, by
$\text{corr}(\bm{c}_{\ell},\bm{c}_{\ell^{\prime}})$, where corr is defined in
Eq. (7). We then average $\text{corr}(\bm{c}_{\ell},\bm{c}_{\ell^{\prime}})$
over all the $K$ matched pairs of states in the two sessions and call it the
centroid similarity. We select the matching that maximizes the centroid
similarity among the $K!$ matchings.

For $K=9$ and $K=10$, we cannot assess all possible pairwise matchings due to
combinatorial explosion. Therefore, we use a greedy search to find an
approximately optimal matching. First, we find the pair of the $\ell$th state
in the first session and the $\ell^{\prime}$th state in the second session
that maximizes $\text{corr}(\bm{c}_{\ell},\bm{c}_{\ell^{\prime}})$. Second, we
select one state from the remaining $K-1$ states in the first session and one
state from the remaining $K-1$ states in the second session such that the
correlation between the two centroids is the largest, and we pair them. We
repeat this procedure until all the $K$ states are matched between the two
sessions.

Once we have determined the final matching between the $K$ states in the first
session and those in the second session, we use the centroid dissimilarity,
defined as $1-(\text{centroid similarity})$, as a measure of discrepancy
between the set of $K$ states in the two sessions. The centroid dissimilarity
ranges between $0$ and $2$. It is equal to $0$ if and only if the set of the
$L$ centroid positions is exactly parallel between the two sessions.

The centroid similarity, $\text{corr}(\bm{c}_{\ell},\bm{c}_{\ell^{\prime}})$,
only compares the direction of the two centroids, $\bm{c}_{\ell}$ and
$\bm{c}_{\ell^{\prime}}$, from the origin. Therefore, we also measured the
discrepancy between the set of $K$ states in the two sessions based on the
Euclidean distance between $\bm{c}_{\ell}$ and $\bm{c}_{\ell^{\prime}}$, given
by

$d(\bm{c}_{\ell},\bm{c}_{\ell^{\prime}})=\left\|\bm{c}_{\ell}-\bm{c}_{\ell^{\prime}}\right\|^{2}.$
(14)

In the verification analysis, we searched for the best matching of the $K$
states between the two sessions by minimizing the average of
$d(\bm{c}_{\ell},\bm{c}_{\ell^{\prime}})$ over the $K$ matched pairs of states
instead of maximizing the average of
$\text{corr}(\bm{c}_{\ell},\bm{c}_{\ell^{\prime}})$. Similar to the case of
using $\text{corr}(\bm{c}_{\ell},\bm{c}_{\ell^{\prime}})$, we did so by the
exhaustive search when $K\in\\{2,\ldots,8\\}$ and by the greedy algorithm when
$K\in\\{9,10\\}$. The dissimilarity obtained using the average of $d$ is equal
to $0$ if and only if the set of the $L$ centroid positions is the same
between the two sessions, and its large value implies a large discrepancy
between the two sessions in terms of the centroid position.

For the coverage time, frequency of appearance, and average lifespan of
states, we compute the total variation (TV) to quantify the difference in the
state-transition dynamics between two sessions. Let $Q_{i}(\ell)$ be the
coverage time, frequency of appearance, or average lifespan for the $\ell$th
state in session $i$. For the notational convenience, we assume without loss
of generality that we have matched the $\ell$th state in session $i$ with the
$\ell$th state in session $j$. For the coverage time of the $\ell$th
microstate, we use the normalized coverage time defined in section 2.6.1 as
$Q_{i}(\ell)$. The TV is defined by

$\delta(Q_{i},Q_{j})=\max_{\ell\in\\{1,2,\ldots,K\\}}\left|Q_{i}(\ell)-Q_{j}(\ell)\right|,$
(15)

where $Q_{i}=\\{Q_{i}(1),\ldots,Q_{i}(K)\\}$.

To quantify the difference between the transition probability matrices for two
sessions $i$ and $j$, denoted by
$P^{(i)}=\left(p^{(i)}_{\ell\ell^{\prime}}\right)$ and
$P^{(j)}=\left(p^{(j)}_{\ell\ell^{\prime}}\right)$, respectively, where
$\ell,\ell^{\prime}\in\\{1,\ldots,K\\}$, we calculate the Frobenius distance
given by

$\left\|P^{(i)}-P^{(j)}\right\|_{F}=\sqrt{\sum_{\ell=1}^{K}\sum_{\ell^{\prime}=1}^{K}\left|P^{(i)}_{\ell\ell^{\prime}}-P^{(j)}_{\ell\ell^{\prime}}\right|^{2}}.$
(16)

####  2.6.3 Permutation test

We hypothesize that the state-transition dynamics estimated from fMRI data is
more consistent between different sessions of the same participant than
between different participants. To test this hypothesis, we compare the
dissimilarity between two sessions originating from the same participant and
the dissimilarity between two sessions originating from different
participants. If the former is smaller than the latter, then the state-
transition dynamics is more reproducible within a participant than between
different participants, supporting the potential ability of state-transition
dynamics to be used for individual fingerprinting.

We measure the dissimilarity between a given pair of sessions in terms of one
of the five observables (i.e., centroid position, distribution of the coverage
time, normalized frequency of appearance of states, distribution of the
average lifespan, or the transition probability matrix). For each observable,
we compare the within-participant dissimilarity and between-participant
dissimilarity using the normalized distance ND combined with the permutation
test [79, 80], which we adapt here for our purpose. Denote by $q(p,s)$ one of
the five observables for participant $p\in\\{1,\ldots,N_{\rm p}\\}$ and
session $s\in\\{1,\ldots,N_{\rm s}\\}$, where $N_{\rm p}=8$ is the number of
participants, and $N_{\rm s}=10$ is the number of sessions per participant. We
define the ND by

$\displaystyle\text{ND}(q)$
$\displaystyle=\frac{\frac{2}{N_{\text{p}}(N_{\text{p}}-1)N_{\text{s}}}\sum\limits_{\begin{subarray}{c}s=1\end{subarray}}^{N_{\text{s}}}\sum\limits_{\begin{subarray}{c}p=1\end{subarray}}^{N_{\text{p}}}\sum\limits_{\begin{subarray}{c}p^{\prime}=1\end{subarray}}^{p-1}\tilde{d}\left(q(p,s),q(p^{\prime},s)\right)}{\frac{2}{N_{\text{p}}N_{\text{s}}(N_{\text{s}}-1)}\sum\limits_{\begin{subarray}{c}p=1\end{subarray}}^{N_{\text{p}}}\sum\limits_{\begin{subarray}{c}s=1\end{subarray}}^{N_{\text{s}}}\sum\limits_{\begin{subarray}{c}s^{\prime}=1\end{subarray}}^{s-1}\tilde{d}\left(q(p,s),q(p,s^{\prime})\right)}$
(17)
$\displaystyle=\frac{(N_{\text{s}}-1)\sum\limits_{\begin{subarray}{c}s=1\end{subarray}}^{N_{\text{s}}}\sum\limits_{\begin{subarray}{c}p=1\end{subarray}}^{N_{\text{p}}}\sum\limits_{\begin{subarray}{c}p^{\prime}=1\end{subarray}}^{p-1}\tilde{d}\left(q(p,s),q(p^{\prime},s)\right)}{(N_{\text{p}}-1)\sum\limits_{\begin{subarray}{c}p=1\end{subarray}}^{N_{\text{p}}}\sum\limits_{\begin{subarray}{c}s=1\end{subarray}}^{N_{\text{s}}}\sum\limits_{\begin{subarray}{c}s^{\prime}=1\end{subarray}}^{s-1}\tilde{d}\left(q(p,s),q(p,s^{\prime})\right)},$

where $\tilde{d}$ denotes the dissimilarity (i.e., the Euclidean distance, TV,
or Frobenius norm, depending on the observable; see section 2.6.2) between two
sessions. The prefactor on the right-hand side on the first line of Eq. (17)
accounts for the normalization; there are
$\frac{N_{\text{p}}(N_{\text{p}}-1)N_{\text{s}}}{2}$ and
$\frac{N_{\text{p}}N_{\text{s}}(N_{\text{s}}-1)}{2}$ terms in the summation on
the numerator and denominator, respectively. Therefore, the numerator of the
right-hand side on the first line of Eq. (17) represents the average
dissimilarity between two sessions obtained from different participants. The
denominator represents the average dissimilarity between two sessions obtained
from the same participant. If the state-transition dynamics are more
consistent among different sessions within the same participant than among
different sessions of different participants, we expect that $\text{ND}(q)>1$.

To statistically test the $\text{ND}(q)$ value, we ran a permutation test
[81]. Specifically, we carried out the following steps.

  * (i)

Shuffle the values of $q$ across all participants and sessions uniformly at
random. This process is equivalent to applying a random permutation on
$\\{q(1,1),q(1,2),\ldots,q(N_{\text{p}},N_{\text{s}})\\}$. We denote the $q$
value for the $s$th session for $p$th participant after the random permutation
by $q^{\prime}(p,s)$. Note that the $q^{\prime}$ value originates from any of
the $N_{\text{p}}$ participants with probability $1/N_{\text{p}}$ and any of
the $N_{\text{s}}$ sessions with probability $1/N_{\text{s}}$.

  * (ii)

Calculate $\text{ND}(q^{\prime})$.

  * (iii)

Repeat steps (i) and (ii) $R$ times. We set $R=10^{4}$.

  * (iv)

The permutation $p$-value is equal to the fraction of the runs among the $R$
runs in which the $\text{ND}(q^{\prime})$ value is larger than the empirical
$\text{ND}(q)$ value.

##  3 Results

###  3.1 Choice of the global signal reduction and clustering methods

We ran the seven clustering methods for each number of clusters,
$K\in\\{2,\ldots,10\\}$, each of the ten sessions, each of the eight
participants, and each of the three global signal removal methods for the DMN
extracted from the MSC data. Then, we calculated the total GEV, i.e.,
$\text{GEV}_{\text{total}}$, for each combination of these variables as a
measure of the quality of clustering. We show the $\text{GEV}_{\text{total}}$
values averaged over all the participants and sessions in Fig. 1(a)–(c), for
each combination of these variations. Each panel of Fig. 1 corresponds to a
treatment of the global signal removal. In all the cases,
$\text{GEV}_{\text{total}}$ increases as $K$ increases. We also find that
$\text{GEV}_{\text{total}}$ is notably larger with the first and second
treatments than the third treatment for all the seven clustering methods and
that $\text{GEV}_{\text{total}}$ is slightly larger under the second than the
first treatment for all values of $K$ and for all the clustering methods.
Because the second treatment of the global signal removal shows the best
performance in terms of clustering quality (i.e., providing the largest
$\text{GEV}_{\text{total}}$), we use the second treatment in the following
analyses.

Figure 1: Performance of estimating discrete states from the DMN extracted
from the MSC data. We show the results for the three treatments of global
signal removal, seven clustering methods, and $K\in\\{2,\ldots,10\\}$.
(a)–(c): Total GEV. (d)–(f): WCSS. (a) and (d): First treatment of the global
signal removal. (b) and (e): Second treatment. (c) and (f): Third treatment.
Each $\text{GEV}_{\text{total}}$ and WCSS value shown is the average over the
eight participants and ten sessions per participant.

We select the three clustering methods with the largest
$\text{GEV}_{\text{total}}$, which are the K-means, TAAHC, and bisecting
K-means. For these three clustering methods, $\text{GEV}_{\text{total}}$ is
around 70% with $K=4$ (K-means: $70.22\pm 2.76\%$ (average $\pm$ standard
deviation calculated on the basis of the 80 sessions of the MSC data), TAAHC:
$68.56\pm 2.98\%$, bisecting K-means: $69.48\pm 2.93\%$) and more than 75%
with $K=7$ (K-means: $77.49\pm 2.07\%$, TAAHC: $76.16\pm 2.29\%$, bisecting
K-means: $75.85\pm 2.31\%$). As references, previous microstate analyses on
EEG data found the $\text{GEV}_{\text{total}}$ values of $70.92\pm 3.65\%$ [9]
and $65.80\pm 4.90\%$ [7] using the K-means clustering, and $69.93\pm 3.58\%$
using TAAHC [9], all with $K=4$. Furthermore, $\text{GEV}_{\text{total}}$
values of $65.03\pm 6.13\%$ and $60.99\pm 5.62\%$ with $K=5$ were reported for
EEG data recorded under eyes-closed and eyes-open conditions, respectively
[65]. A MEG study reported a $\text{GEV}_{\text{total}}$ of $63.97\pm 0.64\%$
using the K-means clustering with $K=10$ [13]. Our present data analysis with
the fMRI data has yielded somewhat larger $\text{GEV}_{\text{total}}$ values
than these studies.

The GEV is based on the similarity in the direction of the $N$-dimensional
fMRI signal, $\bm{x}_{t}$, and the centroid of the cluster, $\bm{c}_{L_{t}}$,
where we remind that $L_{t}$ is the index of the cluster to which $\bm{x}_{t}$
belongs. Therefore, the GEV can be large even if $\bm{x}_{t}$ and
$\bm{c}_{L_{t}}$ are not close to each other. Therefore, we also computed the
WCSS, which is the sum of the distance between $\bm{x}_{t}$ and
$\bm{c}_{L_{t}}$ over all the volumes. We confirmed that the dependence of the
WCSS on the global signal removal method, clustering method, and $K$ is
similar to that with $\text{GEV}_{\text{total}}$ (see Fig. 1(d)–(f)). Note
that a large $\text{GEV}_{\text{total}}$ value implies a good clustering
result, whereas a small WCSS value implies a good clustering result. In
particular, with the WCSS, the second treatment of the global signal removal
is the best among the three treatments, and the best three clustering methods
remain the same, while the GMM performs as equally well as the TAAHC and the
bisecting K-means for the second treatment of the global signal removal (see
Fig. 1(e)).

Therefore, in the remainder of this paper, we further focus our analysis only
on the K-means, TAAHC, and bisecting K-means clustering methods.

###  3.2 Test-retest reliability of the observables of the state-transition
dynamics

We calculated the five observables, i.e., centroid of the clusters, coverage
time, frequency, average lifespan, and transition probability matrix, of the
estimated state-transition dynamics for each of the three selected clustering
methods, each $K$ value, session, and participant. Then, we calculated the
discrepancy in each observable between two sessions. To compare the state-
transition dynamics of different sessions within the same participant, which
we call the within-participant comparison, we calculated the discrepancy in
terms of each observable for each pair of sessions for each participant.
Because there are ten sessions for each of the eight participants, there are
${10\choose 2}\times 8=360$ within-participant comparisons, where $\choose$
represents the binomial coefficient. To compare the state-transition dynamics
between different participants, which we call the between-participant
comparison, we calculated the discrepancy in terms of each observable between
each pair of sessions obtained from different participants. There are
${8\choose 2}\times 10=280$ between-participant comparisons.

We show the distribution of the discrepancy measure for each observable with
$K=4$, separately for the within-participant and between-participant
comparisons, in Figs. 2. Figures 2(a), 2(b), and 2(c) show the results for the
K-means, TAAHC, and bisecting K-means, respectively. We find that the state-
transition dynamics are visually more similar in the within-participant than
between-participant comparisons across all the indices and for all the three
clustering methods when we compare the minimum, maximum, median, first
quartile, and third quartile values of each distribution. For all the three
clustering methods, the gap between the within-participant and between-
participant comparison is apparently the largest for the centroid position
among the five observables. The gap between the within-participant and
between-participant comparisons often looks subtle, in particular for the
coverage time. The results with $K=7$ and $K=10$ are qualitatively the same as
those with $K=4$ (see section S1).

Figure 2: Within-participant and between-participant reproducibility of the
state-transition dynamics with $K=4$ states. (a) K-means. (b) TAAHC. (c)
Bisecting K-means. “Within” and “Between” indicate the within-participant and
between-participant comparisons, respectively. Each box plot shows the
minimum, maximum, median, first quartile, and third quartile of the
measurements. Each dot represents a session. “Centroid” abbreviates the
centroid position, and “Transition prob.” abbreviates the transition
probability matrix.

To test the significance of the difference between the within-participant and
between-participant session-to-session reproducibility of the state-transition
dynamics, we carried out the permutation test. We computed the ND value for
each clustering method, value of $K\in\\{2,\ldots,10\\}$, and observable.
Furthermore, we computed the ND values for $10^{4}$ randomized session-to-
session comparisons. The permutation test concerns whether the ND value for
the original session-to-session comparisons is significantly different from
the ND values for the comparisons between the randomized pairs of sessions.

Figure 3: Distribution of the ND values for the original and randomized
session-to-session comparisons with $K=4$. (a) K-means. (b) TAAHC. (c)
Bisecting K-means. The $p$ values shown are the uncorrected values.

With $K=4$, we show the ND value for the original session-to-session
comparisons and the distribution of the ND values for the randomized sessions
in Fig. 3. Each panel of Fig. 3 shows a combination of the clustering method
and the observable. The vertical dashed lines represent the ND values for the
original session-to-session comparisons. We find that the result of the
permutation test is significant in many cases even after correcting for
multiple comparisons over the three clustering methods, nine values of $K$,
and five observables; see the uncorrected $p$ values in the figure; an
uncorrected $p=0.00037$ corresponds to a Bonferroni corrected $p=0.05$. A
small $p$ value implies that the within-participant session-to-session
reproducibility is higher than the between-participant session-to-session
reproducibility, suggesting the possibility of using the observable for
fingerprinting individuals.

We tabulate the $p$ values from the permutation test for the three clustering
methods, $K\in\\{2,\ldots,10\\}$, and the five observables in Table 2. In the
table, $p=0$ indicates that the ND value for the original session-to-session
comparisons is farther from $1$ than all the $10^{4}$ randomized comparisons,
corresponding to $p<0.5\times 10^{-4}$. The table shows that a majority of the
$p$ values (i.e., 126 out of 135; 93.33%) are smaller than 0.05 (shown with
*). One hundred and seventeen of them (i.e., 86.67% of the 135 comparisons)
remain significant after the Bonferroni correction (shown with ***; equivalent
to $p<0.00037$, uncorrected). Because there are 135 comparisons in the table
and using the Bonferroni correction may be too stringent, we also counted the
cases in which the uncorrected $p$ value is less than $0.001$, shown with **;
there are 119 out of the 135 comparisons (i.e., 88.15%) with $p<0.001$. We
find that the number of significant $p$ values with the K-means and bisecting
K-means is somewhat larger than with the TAAHC (with $p<0.05$, K-means, TAAHC,
and bisecting K-means have 44, 39, and 43 significant comparisons,
respectively). We also find that the $p$ values considerably depend on the
observables. The permutation test result is strongly significant (i.e.,
$p<0.5\times 10^{-4}$) for all the clustering methods and $K$ values for the
centroid position. In contrast, the number of significant combinations of the
clustering method and $K$ value is smallest for the coverage. Lastly, we do
not observe a notable dependence of the permutation test result on $K$.

###  3.3 Robustness tests

For validation, we also estimated the state-transition dynamics for the whole-
brain network extracted from the MSC data. We show the results for the
permutation test in Table 3. The results are similar to those for the DMN. In
particular, the centroid position is the most effective among the five
observables at distinguishing between the within-participant and between-
participant comparisons, and the coverage is the least effective.

We also ran the same permutation test for a whole-brain network obtained from
the resting-state HCP data. The results, shown in Table 4, are similar to
those for the DMN and whole-brain networks extracted from the MSC data.
However, the permutation test results were stronger for the HCP than MSC data
(i.e., a larger number of significant comparisons among the 135 comparisons).
In particular, the frequency, lifespan, and the transition probability matrix
as well as the centroid position yielded the smallest possible $p$ value
(i.e., $p<0.5\times 10^{-4}$) for all pairs of the clustering method and $K$
value.

Table 2: Results of the permutation test for the DMN extracted from the MSC data. The entries with 0 mean that $p<0.5\cdot 10^{-4}$, uncorrected. *: $p<0.05$, uncorrected; **: $p<0.001$, uncorrected; ***: $p<0.05$, Bonferroni corrected (which is equivalent to $p<0.00037$, uncorrected). We remark that “Centroid” and ”Transition prob.” abbreviate the centroid’s position and the transition probability matrix, respectively.  |  $K$ |  Centroid |  Coverage |  Frequency |  Lifespan |  Transition prob.
---|---|---|---|---|---|---
K-means |  2 |  $0^{***}$ |  0.1975 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
TAAHC |  2 |  $0^{***}$ |  $0.1077$ |  $0.0001^{***}$ |  $0^{***}$ |  $0.0001^{***}$
|  3 |  $0^{***}$ |  $0.0001^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0.3659$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0.0928$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.0162^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  0.5488 |  $0^{***}$ |  $0.0003^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0.0642$ |  $0^{**}$ |  $0.0001^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0.0942$ |  $0.0001^{***}$ |  $0^{**}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.0132^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
Bisecting K-means |  2 |  $0^{***}$ |  0.0568 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  0.0628 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0.0114^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.0102^{*}$ |  $0^{***}$ |  $0.0003^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.0003^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0.0017^{*}$ |  $0.0004^{**}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0.0112^{*}$ |  $0.0009^{**}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.0162^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
Table 3: Results of the permutation test for the whole-brain network extracted from the MSC data. The entries with 0 mean that $p<0.5\cdot 10^{-4}$, uncorrected. *: $p<0.05$, uncorrected; **: $p<0.001$, uncorrected; ***: $p<0.05$, Bonferroni corrected. |  $K$ |  Centroid |  Coverage |  Frequency |  Lifespan |  Transition prob.
---|---|---|---|---|---|---
K-means |  2 |  $0^{***}$ |  $0.045^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0.0193^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0.0414^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.0002^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.0008^{**}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
TAAHC |  2 |  $0^{***}$ |  $0.1747$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0.0167^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  0.4553 |  $0^{***}$ |  $0.0001^{***}$ |  $0^{**}$
|  5 |  $0^{***}$ |  $0.0288^{*}$ |  $0.0001^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.5063$ |  $0.0002^{***}$ |  $0^{**}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.7551$ |  $0.0008^{**}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  0.1887 |  $0.0001^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0.0874$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.0426^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
Bisecting K-means |  2 |  $0^{***}$ |  $0.0283^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0.0597$ |  $0^{***}$ |  $0.0007^{**}$ |  $0.0012^{*}$
|  4 |  $0^{***}$ |  $0.0402^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  0.2698 |  $0.0001^{***}$ |  $0.0003^{***}$ |  $0.065$
|  6 |  $0^{***}$ |  $0.0341^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.5185$ |  $0.0002^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  0.0886 |  $0.0027^{*}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.0012^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
Table 4: Results of the permutation test for the whole-brain network extracted from the HCP data. The entries with 0 mean that the $p<0.5\cdot 10^{-4}$, uncorrected. *: $p<0.05$, uncorrected; **: $p<0.001$, uncorrected; ***: $p<0.05$, Bonferroni corrected. |  $K$ |  Centroid |  Coverage |  Frequency |  Lifespan |  Transition prob.
---|---|---|---|---|---|---
K-means |  2 |  $0^{***}$ |  $0.0308^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
TAAHC |  2 |  $0^{***}$ |  $0.0001^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0.0016^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0.0305^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0.0297^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.0164^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.2864$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0.0699$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0.0203^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.0764$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
Bisecting K-means |  2 |  $0^{***}$ |  $0.0015^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0.0531$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0.0078^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.0009^{**}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.005^{**}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$

Lastly, as we noted earlier, our main definition of the centroid dissimilarity
relies on the (dis)similarity between $\bm{c}_{L_{t}}$ and $\bm{x}_{t}$ only
in terms of the direction. Therefore, we reran the permutation test by
replacing the centroid (dis)similarity by the WCSS to measure the average
distance between $\bm{c}_{L_{t}}$ and $\bm{x}_{t}$. This change not only
affected the discrepancy measure between two sessions in terms of the centroid
position but also the discrepancy between pairs of sessions in terms of the
other four observables (i.e., coverage time, frequency of appearance of each
state, average lifespan, and transition probability matrix). This is because
changing the discrepancy measure for cluster centroids affects how the set of
centroids (and therefore clusters) is matched between two given sessions. We
confirmed that the permutation test results with the WCSS is similar to those
with $\text{GEV}_{\text{total}}$ (see section S2). In particular, the $p$
values were overall small, the results tended to be more significant for the
K-means and bisecting K-means than for the TAAHC, and for the centroid
position and the transition probability matrix than for the other three
observables.

##  4 Discussion

We carried out a comparative study of methods to cluster volumes of the fMRI
to extract time series of the system’s state, akin to microstate analysis for
EEG and MEG data, for each recording session. We found that aggregating the
symmetrically located ROIs into one ROI and then conducting the global signal
removal yielded a high accuracy of clustering in terms of the total GEV and
WCSS. We obtained total GEV values that are somewhat larger than those
obtained in previous studies for EEG microstate analysis [9, 7, 65, 13], which
suggests that fMRI state-transition dynamics analysis may be promising.
Furthermore, by carrying over the three clustering methods yielding the best
clustering performance to a test-retest reliability analysis, we found that,
for different fMRI data sets and different networks, test-retest reliability
was higher in the within-participant comparison than the between-participant
comparison. This result held true for most combinations of the number of
clusters, $K\in\\{2,\ldots,10\\}$, and index quantifying the estimated state-
transition dynamics. We also found that the K-means clustering yielded the
highest test-retest reliability among the three clustering methods. The
present results suggest that clustering-based analysis of state-transition
dynamics, which is substantially simpler than the hidden Markov model [14, 23,
24, 25, 21, 82, 27, 22, 12, 13] and the energy landscape analysis [28, 29, 30,
31], may be a sufficiently competitive method to derive state-transition
dynamics in fMRI data.

The microstate analysis was originally proposed for EEG data [83, 8, 5, 84, 6,
7, 11]. Microstates in EEG data are typically of the order of 100 ms. One
cannot directly associate the discrete states estimated from fMRI data with
EEG or MEG microstates because the time resolution of fMRI data is much lower
than 100 ms; a typical TR is approximately between 1 to 3 seconds.
Furthermore, the typical duration of a discrete state is longer than one TR.
For example, the average lifespan of a state was $3.3$ TR, $2.5$ TR, and $2.2$
TR when we estimated four, seven, and ten states, respectively, for the DMN
extracted from the MSC data. Therefore, cognitively or physiologically
relevant discrete states estimated for fMRI data [19, 66, 24] may be different
from those captured by microstates in EEG and MEG data. However, promising
correspondences between EEG microstates and fMRI states have been reported
[15, 85, 19, 86]. Analyzing simultaneously recorded EEG-fMRI data may further
reveal connection between EEG microstates and discrete states for fMRI data
[15, 87, 88, 89, 90, 18, 91, 92].

We examined test-retest reliability of discrete states estimated by clustering
activity pattern vectors of fMRI data. In contrast, various previous studies
estimated discrete states by clustering functional networks from fMRI data
[36, 37, 38, 20, 14, 19, 39, 23, 25, 21]. Our methods of test-retest
reliability analysis do not depend on how the discrete states are estimated
and therefore are applicable to the case of state-transition dynamics of
functional networks. To the best of our knowledge, no work has systematically
compared the reliability between state-transition dynamics estimated for
spatial activity patterns or their vectorized versions and those estimated for
functional networks from the same fMRI data. Such a comparative analysis may
better inform us either activity patterns or functional networks are more
powerful biomarkers than the other when combined with state-transition
dynamics modeling. In a similar vein, the aforementioned studies pursuing
similarity between EEG microstates and fMRI dynamic states have been confined
to the case in which fMRI dynamic states are estimated from dynamics of
functional connectivity, not dynamics of activity patterns. These topics
warrant future work.

In EEG microstate analysis, it is common to generate global microstate maps,
which is to determine a given number of microstates by clustering candidate
EEG maps obtained from different participants altogether. Then, one matches
the obtained microstate maps shared by all the participants to the individual
EEG maps from the individual participants to determine the microstate dynamics
for each participant [93, 15, 9, 66, 65]. For EEG data, this approach has been
shown to accrue higher reliability than microstate maps estimated separately
for individual participants [9]. Nevertheless, we have estimated the states
separately for individual participants (and for individual sessions) in the
present study. This is because, for fMRI data, one often estimates state
dynamics separately for each individual, which allows one to study subject
variability of the estimated state dynamics or to exploit it [19, 94, 25, 21].
In contrast, pooling fMRI data from different participants to generate across-
participant templates of discrete states is also a common practice [32, 94,
19, 66, 35]. In fact, one can run our test-retest reliability analysis even if
we estimate the templates of the discrete states shared by all participants,
with the exception of the cluster centroid, $\bm{c}_{\ell}$, as an observable
of the estimated state dynamics; if we use a shared template, $\bm{c}_{\ell}$
is the same for all sessions and individuals and therefore one cannot compare
its reliability within versus between participants. We point out that
comparison of the reliability between shared templates and individualized
templates of discrete states for fMRI data, as was done for EEG data [9], is
underexplored.

We ran a permutation test to statistically compare the within-participant and
between-participant test-retest reliability. This permutation test is an
adaptation of what we recently developed for energy landscape analysis [50] to
the case of clustering-based state-transition dynamics. This method is not
limited to fMRI data. It is straightforward to use it for EEG and MEG
microstate data analysis obtained from multiple participants and multiple
sessions per participant. Our code is publicly available on Github [49]. The
only requirement is to define observables and to be able to measure the
discrepancy in the observable between an arbitrary pair of sessions. Assessing
test-retest reliability in EEG [9, 6, 10, 11] and MEG [95, 96] data using this
technique as well as furthering the application to fMRI data in health and
disease may be fruitful.

## Acknowledgements

T.W. acknowledges support from the Japan Society for Promotion of Sciences
(TW, 19H03535, 21H05679, 23H04217) N.M. acknowledges support from the Japan
Science and Technology Agency (JST) Moonshot R&D (under grant no. JPMJMS2021),
the National Science Foundation (under grant no. 2204936), and JSPS KAKENHI
(under grant nos. JP 21H04595 and 23H03414). Two publicly available data sets
were used in this work. The first data set was provided by the Midnight Scan
Club (MSC) project, funded by NIH Grants NS088590, TR000448 (NUFD), MH104592
(DJG), and HD087011 (to the Intellectual and Developmental Disabilities
Research Center at Washington University); the Jacobs Foundation (NUFD); the
Child Neurology Foundation (NUFD); the McDonnell Center for Systems
Neuroscience (NUFD, BLS); the Mallinckrodt Institute of Radiology (NUFD); the
Hope Center for Neurological Disorders (NUFD, BLS, SEP); and Dart Neuroscience
LLC. This data was obtained from the OpenfMRI database. Its accession number
is ds000224. The second data set was provided by the Human Connectome Project,
WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil
Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that
support the NIH Blueprint for Neuroscience Research; and by the McDonnell
Center for Systems Neuroscience at Washington University. The authors also
acknowledge support provided by the Center for Computational Research at the
University at Buffalo for the processing of HCP data.

## References

  * [1] M. E. Raichle, A. M. MacLeod, A. Z. Snyder, W. J. Powers, D. A. Gusnard, G. L. Shulman, A default mode of brain function, Proceedings of the National Academy of Sciences of the United States of America 98 (2001) 676–682.
  * [2] M. I. Rabinovich, P. Varona, A. I. Selverston, H. D. Abarbanel, Dynamical principles in neuroscience, Reviews of Modern Physics 78 (2006) 1213.
  * [3] G. Deco, G. Tononi, M. Boly, M. L. Kringelbach, Rethinking segregation and integration: contributions of whole-brain modelling, Nature Reviews Neuroscience 16 (2015) 430–439.
  * [4] A. Avena-Koenigsberger, B. Misic, O. Sporns, Communication dynamics in complex brain networks, Nature Reviews Neuroscience 19 (2018) 17–33.
  * [5] A. Khanna, A. Pascual-Leone, C. M. Michel, F. Farzan, Microstates in resting-state EEG: Current status and future directions, Neuroscience & Biobehavioral Reviews 49 (2015) 105–113.
  * [6] C. M. Michel, T. Koenig, EEG microstates as a tool for studying the temporal dynamics of whole-brain neuronal networks: A review, NeuroImage 180 (2018) 577–593.
  * [7] F. Von Wegner, P. Knaut, H. Laufs, EEG microstate sequences from different clustering algorithms are information-theoretically invariant, Frontiers in Computational Neuroscience 12 (2018) 70.
  * [8] T. Koenig, L. Prichep, D. Lehmann, P. V. Sosa, E. Braeker, H. Kleinlogel, R. Isenhart, E. R. John, Millisecond by millisecond, year by year: Normative EEG microstates and developmental stages, NeuroImage 16 (2002) 41–48.
  * [9] A. Khanna, A. Pascual-Leone, F. Farzan, Reliability of resting-state microstate features in electroencephalography, PLoS ONE 9 (2014) e114163.
  * [10] J. Liu, J. Xu, G. Zou, Y. He, Q. Zou, J.-H. Gao, Reliability and individual specificity of EEG microstate characteristics, Brain Topography 33 (2020) 438–449.
  * [11] K. Zhang, W. Shi, C. Wang, Y. Li, Z. Liu, T. Liu, J. Li, X. Yan, Q. Wang, Z. Cao, et al., Reliability of EEG microstate analysis at different electrode densities during propofol-induced transitions of brain states, NeuroImage 231 (2021) 117861.
  * [12] L. Tait, J. Zhang, +microstate: A MATLAB toolbox for brain microstate analysis in sensor and cortical EEG/MEG, NeuroImage 258 (2022) 119346.
  * [13] L. Tait, J. Zhang, MEG cortical microstates: spatiotemporal characteristics, dynamic functional connectivity and stimulus-evoked responses, NeuroImage 251 (2022) 119006.
  * [14] A. P. Baker, M. J. Brookes, I. A. Rezek, S. M. Smith, T. Behrens, P. J. Probert Smith, M. Woolrich, Fast transient networks in spontaneous human brain activity, eLife 3 (2014) e01867.
  * [15] J. Britz, D. Van De Ville, C. M. Michel, Bold correlates of EEG topography reveal rapid resting-state network dynamics, NeuroImage 52 (2010) 1162–1170.
  * [16] S. Schwab, T. Koenig, Y. Morishima, T. Dierks, A. Federspiel, K. Jann, Discovering frequency sensitive thalamic nuclei from EEG microstate informed resting state fMRI, NeuroImage 118 (2015) 368–375.
  * [17] M. Case, H. Zhang, J. Mundahl, Y. Datta, S. Nelson, K. Gupta, B. He, Characterization of functional brain activity and connectivity using EEG and fMRI in patients with sickle cell disease, NeuroImage: Clinical 14 (2017) 1–17.
  * [18] L. Bréchet, D. Brunet, G. Birot, R. Gruetter, C. M. Michel, J. Jorge, Capturing the spatiotemporal dynamics of self-generated, task-initiated thoughts with EEG and fMRI, NeuroImage 194 (2019) 82–92.
  * [19] E. A. Allen, E. Damaraju, S. M. Plis, E. B. Erhardt, T. Eichele, V. D. Calhoun, Tracking whole-brain connectivity dynamics in the resting state, Cerebral Cortex 24 (2014) 663–676.
  * [20] V. D. Calhoun, R. Miller, G. Pearlson, T. Adalı, The chronnectome: Time-varying connectivity networks as the next frontier in fMRI data discovery, Neuron 84 (2014) 262–274.
  * [21] S. F. Nielsen, M. N. Schmidt, K. H. Madsen, M. Mørup, Predictive assessment of models for dynamic functional connectivity, NeuroImage 171 (2018) 116–134.
  * [22] T. Ezaki, Y. Himeno, T. Watanabe, N. Masuda, Modelling state-transition dynamics in resting-state brain signals by the hidden Markov and Gaussian mixture models, European Journal of Neuroscience 54 (2021) 5404–5416.
  * [23] S. Ryali, K. Supekar, T. Chen, J. Kochalka, W. Cai, J. Nicholas, A. Padmanabhan, V. Menon, Temporal dynamics and developmental maturation of salience, default and central-executive network interactions revealed by variational Bayes hidden Markov modeling, PLoS Computational Biology 12 (2016) e1005138.
  * [24] D. Vidaurre, S. M. Smith, M. W. Woolrich, Brain network dynamics are hierarchically organized in time, Proceedings of the National Academy of Sciences of the United States of America 114 (2017) 12827–12832.
  * [25] J. Taghia, S. Ryali, T. Chen, K. Supekar, W. Cai, V. Menon, Bayesian switching factor analysis for estimating time-varying functional connectivity in fMRI, NeuroImage 155 (2017) 271–290.
  * [26] M. J. Brookes, M. J. Groom, L. Liuzzi, R. M. Hill, H. J. Smith, P. M. Briley, E. L. Hall, B. A. Hunt, L. E. Gascoyne, M. J. Taylor, et al., Altered temporal stability in dynamic neural networks underlies connectivity changes in neurodevelopment, NeuroImage 174 (2018) 563–575.
  * [27] D. Vidaurre, A new model for simultaneous dimensionality reduction and time-varying functional connectivity estimation, PLoS Computational Biology 17 (2021) e1008580.
  * [28] T. Watanabe, S. Hirose, H. Wada, Y. Imai, T. Machida, I. Shirouzu, S. Konishi, Y. Miyashita, N. Masuda, Energy landscapes of resting-state brain networks, Frontiers in Neuroinformatics 8 (2014) 12.
  * [29] T. Watanabe, N. Masuda, F. MEGumi, R. Kanai, G. Rees, Energy landscape and dynamics of brain activity during human bistable perception, Nature Communications 5 (2014) 4765.
  * [30] T. Ezaki, T. Watanabe, M. Ohzeki, N. Masuda, Energy landscape analysis of neuroimaging data, Philosophical Transactions of the Royal Society A 375 (2017) 20160287.
  * [31] T. Ezaki, M. Sakaki, T. Watanabe, N. Masuda, Age-related changes in the ease of dynamical transitions in human brain activity, Human Brain Mapping 39 (2018) 2673–2688.
  * [32] X. Liu, C. Chang, J. H. Duyn, Decomposition of spontaneous brain activity into distinct fMRI co-activation patterns, Frontiers in Systems Neuroscience 7 (2013) 101.
  * [33] X. Liu, J. H. Duyn, Time-varying functional network information extracted from brief instances of spontaneous brain activity, Proceedings of the National Academy of Sciences of the United States of America 110 (2013) 4392–4397.
  * [34] X. Liu, N. Zhang, C. Chang, J. H. Duyn, Co-activation patterns in resting-state fMRI signals, NeuroImage 180 (2018) 485–494.
  * [35] J.-J. Paakki, J. S. Rahko, A. Kotila, M.-L. Mattila, H. Miettunen, T. M. Hurtig, K. K. Jussila, S. Kuusikko-Gauffin, I. K. Moilanen, O. Tervonen, et al., Co-activation pattern alterations in autism spectrum disorder–a volume-wise hierarchical clustering fMRI study, Brain and Behavior 11 (2021) e02174.
  * [36] Ü. Sakoğlu, G. D. Pearlson, K. A. Kiehl, Y. M. Wang, A. M. Michael, V. D. Calhoun, A method for evaluating dynamic functional network connectivity and task-modulation: application to schizophrenia, Magnetic Resonance Materials in Physics, Biology and Medicine 23 (2010) 351–366.
  * [37] R. M. Hutchison, T. Womelsdorf, J. S. Gati, S. Everling, R. S. Menon, Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques, Human Brain Mapping 34 (2013) 2154–2177.
  * [38] N. Leonardi, J. Richiardi, M. Gschwind, S. Simioni, J.-M. Annoni, M. Schluep, P. Vuilleumier, D. Van De Ville, Principal components of functional connectivity: A new approach to study dynamic brain connectivity during rest, NeuroImage 83 (2013) 937–950.
  * [39] A. Abrol, C. Chaze, E. Damaraju, V. D. Calhoun, The chronnectome: Evaluating replicability of dynamic connectivity patterns in 7500 resting fMRI datasets, in: 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), IEEE, 2016, pp. 5571–5574.
  * [40] E. S. Finn, X. Shen, D. Scheinost, M. D. Rosenberg, J. Huang, M. M. Chun, X. Papademetris, R. T. Constable, Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity, Nature Neuroscience 18 (2015) 1664–1671.
  * [41] E. M. Gordon, T. O. Laumann, A. W. Gilmore, D. J. Newbold, D. J. Greene, J. J. Berg, M. Ortega, C. Hoyt-Drazen, C. Gratton, H. Sun, et al., Precision functional mapping of individual human brains, Neuron 95 (2017) 791–807.
  * [42] E. Amico, J. Goñi, The quest for identifiability in human functional connectomes, Scientific Reports 8 (2018) 8254.
  * [43] S. Noble, D. Scheinost, R. T. Constable, A decade of test-retest reliability of functional connectivity: A systematic review and meta-analysis, NeuroImage 203 (2019) 116157.
  * [44] M. Venkatesh, J. Jaja, L. Pessoa, Comparing functional connectivity matrices: A geometry-aware approach applied to participant identification, NeuroImage 207 (2020) 116398.
  * [45] B. Chiêm, K. Abbas, E. Amico, D. A. Duong-Tran, F. Crevecoeur, J. Goñi, Improving functional connectome fingerprinting with degree-normalization, Brain Connectivity 12 (2022) 180–192.
  * [46] C. Zhang, S. A. Baum, V. R. Adduru, B. B. Biswal, A. M. Michael, Test-retest reliability of dynamic functional connectivity in resting state fMRI, NeuroImage 183 (2018) 907–918.
  * [47] X. Zhang, J. Liu, Y. Yang, S. Zhao, L. Guo, J. Han, X. Hu, Test–retest reliability of dynamic functional connectivity in naturalistic paradigm functional magnetic resonance imaging, Human Brain Mapping 43 (2022) 1463–1476.
  * [48] Y. Long, X. Ouyang, C. Yan, Z. Wu, X. Huang, W. Pu, H. Cao, Z. Liu, L. Palaniyappan, Evaluating test–retest reliability and sex-/age-related effects on temporal clustering coefficient of dynamic functional brain networks, Human Brain Mapping 44 (2023) 2191–2208.
  * [49] S. Islam, Functional MRI state transition dynamics, github repository, accessed: August 5, 2023 (2023).
URL https://github.com/sislam99/fmri_state_transition_dynamics

  * [50] P. Khanra, J. Nakuci, S. Muldoon, T. Watanabe, N. Masuda, Reliability of energy landscape analysis of resting-state functional MRI data, arXiv preprint arXiv:2305.19573 (2023).
  * [51] D. A. Fair, A. L. Cohen, J. D. Power, N. U. Dosenbach, J. A. Church, F. M. Miezin, B. L. Schlaggar, S. E. Petersen, Functional brain networks develop from a “local to distributed” organization, PLoS Computational Biology 5 (2009) e1000381.
  * [52] 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 72 (2011) 665–678.
  * [53] T. Watanabe, G. Rees, Brain network dynamics in high-functioning individuals with autism, Nature Communications 8 (2017) 1–14.
  * [54] D. C. Van Essen, S. M. Smith, D. M. Barch, T. E. Behrens, E. Yacoub, K. Ugurbil, W.-M. H. Consortium, et al., The WU-Minn Human Connectome Project: An overview, NeuroImage 80 (2013) 62–79.
  * [55] M. F. Glasser, S. N. Sotiropoulos, J. A. Wilson, T. S. Coalson, B. Fischl, J. L. Andersson, J. Xu, S. Jbabdi, M. Webster, J. R. Polimeni, et al., The minimal preprocessing pipelines for the Human Connectome Project, NeuroImage 80 (2013) 105–124.
  * [56] M. Jenkinson, P. Bannister, M. Brady, S. Smith, Improved optimization for the robust and accurate linear registration and motion correction of brain images, NeuroImage 17 (2002) 825–841.
  * [57] J. D. Power, K. A. Barnes, A. Z. Snyder, B. L. Schlaggar, S. E. Petersen, Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion, NeuroImage 59 (2012) 2142–2154.
  * [58] G. C. Burgess, S. Kandala, D. Nolan, T. O. Laumann, J. D. Power, B. Adeyemo, M. P. Harms, S. E. Petersen, D. M. Barch, Evaluation of denoising strategies to address motion-correlated artifacts in resting-state functional magnetic resonance imaging data from the Human Connectome Project, Brain Connectivity 6 (2016) 669–680.
  * [59] J. S. Siegel, A. Mitra, T. O. Laumann, B. A. Seitzman, M. Raichle, M. Corbetta, A. Z. Snyder, Data quality influences observed links between functional connectivity and behavior, Cerebral Cortex 27 (2017) 4492–4502.
  * [60] N. Tzourio-Mazoyer, B. Landeau, D. Papathanassiou, F. Crivello, O. Etard, N. Delcroix, B. Mazoyer, M. Joliot, Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain, NeuroImage 15 (2002) 273–289.
  * [61] A. Schaefer, R. Kong, E. M. Gordon, T. O. Laumann, X.-N. Zuo, A. J. Holmes, S. B. Eickhoff, B. T. Yeo, Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity MRI, Cerebral Cortex 28 (2018) 3095–3114.
  * [62] K. Murphy, M. D. Fox, Towards a consensus regarding global signal regression for resting state functional connectivity MRI, NeuroImage 154 (2017) 169–173.
  * [63] J. Li, T. Bolt, D. Bzdok, J. S. Nomi, B. Yeo, R. N. Spreng, L. Q. Uddin, Topography and behavioral relevance of the global signal in the human brain, Scientific Reports 9 (2019) 14286.
  * [64] K. M. Aquino, B. D. Fulcher, L. Parkes, K. Sabaroedin, A. Fornito, Identifying and removing widespread signal deflections from fMRI data: Rethinking the global signal regression problem, NeuroImage 212 (2020) 116614.
  * [65] A. P. Zanesco, B. G. King, A. C. Skwara, C. D. Saron, Within and between-person correlates of the temporal dynamics of resting EEG microstates, NeuroImage 211 (2020) 116631.
  * [66] A. Abrol, E. Damaraju, R. L. Miller, J. M. Stephen, E. D. Claus, A. R. Mayer, V. D. Calhoun, Replicability of time-varying connectivity patterns in large resting state fMRI samples, NeuroImage 163 (2017) 160–176.
  * [67] M. Khosla, K. Jamison, G. H. Ngo, A. Kuceyeski, M. R. Sabuncu, Machine learning in resting-state fMRI analysis, Magnetic Resonance Imaging 64 (2019) 101–121.
  * [68] D. Arthur, S. Vassilvitskii, K-means++ the advantages of careful seeding, in: Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, 2007, pp. 1027–1035.
  * [69] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, et al., Scikit-learn: Machine learning in python, The Journal of Machine Learning Research 12 (2011) 2825–2830.
  * [70] L. Kaufman, P. J. Rousseeuw, Finding Groups in Data: An Introduction to Cluster Analysis, John Wiley & Sons, Hoboken, NJ, 2009.
  * [71] R. Tibshirani, G. Walther, Cluster validation by prediction strength, Journal of Computational and Graphical Statistics 14 (2005) 511–528.
  * [72] M. M. Murray, D. Brunet, C. M. Michel, Topographic ERP analyses: A step-by-step tutorial review, Brain Topography 20 (2008) 249–264.
  * [73] A. T. Poulsen, A. Pedroni, N. Langer, L. K. Hansen, Microstate EEGlab toolbox: An introductory guide, BioRxiv (2018) 289850.
  * [74] R. D. Pascual-Marqui, C. M. Michel, D. Lehmann, Segmentation of brain electrical activity into microstates: model estimation and validation, IEEE Transactions on Biomedical Engineering 42 (1995) 658–665.
  * [75] M. Steinbach, G. Karypis, V. Kumar, A comparison of document clustering techniques, Tech. Rep. 00-034, Department of Computer Science and Egineering, University of Minnesota, https://hdl.handle.net/11299/215421. Accessed: October 04, 2022 (2000).
  * [76] C. M. Bishop, N. M. Nasrabadi, Pattern Recognition and Machine Learning, Vol. 4, Springer, New York, NY, 2006.
  * [77] B. G. Lindsay, Mixture models: theory, geometry, and applications, in: NSF-CBMS Regional Conference Series in Probability and Statistics, Vol. 5, Institute of Mathematical Statistics, The United States of America, 1995, pp. i–163.
  * [78] A. K. Jain, R. C. Dubes, Algorithms for Clustering Data, Prentice-Hall, Inc., Englewood Cliffs, NJ, 1988.
  * [79] S. Bernardi, M. K. Benna, M. Rigotti, J. Munuera, S. Fusi, C. D. Salzman, The geometry of abstraction in the hippocampus and prefrontal cortex, Cell 183 (2020) 954–967.
  * [80] Y. Liu, S. L. Brincat, E. K. Miller, M. E. Hasselmo, A geometric characterization of population coding in the prefrontal cortex and hippocampus during a paired-associate learning task, Journal of Cognitive Neuroscience 32 (2020) 1455–1465.
  * [81] E. Maris, R. Oostenveld, Nonparametric statistical testing of EEG-and MEG-data, Journal of Neuroscience Methods 164 (2007) 177–190.
  * [82] D. Vidaurre, L. T. Hunt, A. J. Quinn, B. A. Hunt, M. J. Brookes, A. C. Nobre, M. W. Woolrich, Spontaneous cortical activity transiently organises into frequency specific phase-coupling networks, Nature Communications 9 (2018) 2987\.
  * [83] D. Lehmann, H. Ozaki, I. Pal, EEG alpha map series: brain micro-states by space-oriented adaptive segmentation, Electroencephalography and Clinical Neurophysiology 67 (1987) 271–288.
  * [84] F. Hatz, M. Hardmeier, H. Bousleiman, S. Rüegg, C. Schindler, P. Fuhr, Reliability of functional connectivity of electroencephalography applying microstate-segmented versus classical calculation of phase lag index, Brain Connectivity 6 (2016) 461–469.
  * [85] H. Yuan, V. Zotev, R. Phillips, W. C. Drevets, J. Bodurka, Spatiotemporal dynamics of the brain at rest—exploring EEG microstates as electrophysiological signatures of bold resting state networks, NeuroImage 60 (2012) 2062–2072.
  * [86] R. Abreu, J. Jorge, A. Leal, T. Koenig, P. Figueiredo, EEG microstates predict concurrent fMRI dynamic functional connectivity states, Brain Topography 34 (2021) 41–55.
  * [87] C. Chang, Z. Liu, M. C. Chen, X. Liu, J. H. Duyn, EEG correlates of time-varying bold functional connectivity, NeuroImage 72 (2013) 227–236.
  * [88] M. G. Preti, N. Leonardi, F. I. Karahanoğlu, F. Grouiller, M. Genetti, M. Seeck, S. Vulliemoz, D. Van De Ville, Epileptic network activity revealed by dynamic functional connectivity in simultaneous EEG-fMRI, in: 2014 IEEE 11th International Symposium on Biomedical Imaging (ISBI), IEEE, 2014, pp. 9–12.
  * [89] R. F. Ahmad, A. S. Malik, N. Kamel, F. Reza, J. M. Abdullah, Simultaneous EEG-fMRI for working memory of the human brain, Australasian Physical & Engineering Sciences in Medicine 39 (2016) 363–378.
  * [90] T. Keinänen, S. Rytky, V. Korhonen, N. Huotari, J. Nikkinen, O. Tervonen, J. M. Palva, V. Kiviniemi, Fluctuations of the EEG-fMRI correlation reflect intrinsic strength of functional connectivity in default mode network, Journal of Neuroscience Research 96 (2018) 1689–1698.
  * [91] D. J. Lurie, D. Kessler, D. S. Bassett, R. F. Betzel, M. Breakspear, S. Kheilholz, A. Kucyi, R. Liégeois, M. A. Lindquist, A. R. McIntosh, et al., Questions and controversies in the study of time-varying functional connectivity in resting fMRI, Network Neuroscience 4 (2020) 30–69.
  * [92] C. Mulert, Simultaneous EEG and fMRI: towards the characterization of structure and dynamics of brain networks, Dialogues in Clinical Neuroscience 15 (2013) 381–386.
  * [93] D. Lehmann, P. L. Faber, S. Galderisi, W. M. Herrmann, T. Kinoshita, M. Koukkou, A. Mucci, R. D. Pascual-Marqui, N. Saito, J. Wackermann, et al., EEG microstate duration and syntax in acute, medication-naive, first-episode schizophrenia: a multi-center study, Psychiatry Research: NeuroImaging 138 (2005) 141–156.
  * [94] B. Rashid, E. Damaraju, G. D. Pearlson, V. D. Calhoun, Dynamic connectivity states estimated from resting fMRI identify differences among Schizophrenia, bipolar disorder, and healthy control subjects, Frontiers in Human Neuroscience 8 (2014) 897.
  * [95] P. Garcés, M. C. Martín-Buro, F. Maestú, Quantifying the test-retest reliability of magnetoencephalography resting-state functional connectivity, Brain Connectivity 6 (2016) 448–460.
  * [96] F. T. Candelaria-Cook, M. E. Schendel, C. J. Ojeda, J. R. Bustillo, J. M. Stephen, Reduced parietal alpha power and psychotic symptoms: Test-retest reliability of resting-state magnetoencephalography in schizophrenia and healthy controls, Schizophrenia Research 215 (2020) 229–240.

Supplementary Information for
State-transition dynamics of resting-state functional magnetic resonance
imaging data: Model comparison and test-to-retest analysis
Saiful Islam1, Pitambar Khanra2, Johan Nakuci4, Sarah F. Muldoon1,2,3,
Takamitsu Watanabe5, Naoki Masuda1,2∗

1Institute for Artificial Intelligence and Data Science, University at
Buffalo, State University of New York at Buffalo, USA
2Department of Mathematics, University at Buffalo, State University of New
York, USA
3Neuroscience Program, University at Buffalo, State University of New York at
Buffalo, USA
4School of Psychology, Georgia Institute of Technology, USA
5International Research Centre for Neurointelligence, The University of Tokyo
Institutes for Advanced Study, Tokyo, Japan
∗Corresponding author<EMAIL_ADDRESS>
##  S1 Within-participant and between-participant reproducibility of the
state-transition dynamics with $K=7$ and $K=10$

We show the distributions of the discrepancy in terms of the five observables
for the within-participant and between-participant comparisons in Figs. S1 and
S2 for the number of clusters $K=7$ and $K=10$, respectively.

Figure S1: Within-participant and between-participant reproducibility of the
state-transition dynamics with $K=7$ states. (a) K-means. (b) TAAHC. (c)
Bisecting K-means. “Within” and “Between” indicate the within-participant and
between-participant comparisons, respectively. Each box plot shows the
minimum, maximum, median, first quartile, and third quartile of the
measurements. Each dot represents a session. “Centroid” abbreviates the
centroid position, and “Transition prob.” abbreviates the transition
probability matrix. Figure S2: Within-participant and between-participant
reproducibility of the state-transition dynamics with $K=10$ states. (a)
K-means. (b) TAAHC. (c) Bisecting K-means.

##  S2 Robustness test with the Euclidean distance used as dissimilarity
between sessions

We mainly used the dissimilarity based on the
$\text{corr}(\bm{c}_{\ell},\bm{c}_{\ell^{\prime}})$ defined by Eq. (7) in the
main text to match and compare the set of the centroid positions,
$\\{\bm{c}_{1},\ldots,\bm{c}_{L}\\}$, between pairs of sessions. In this
section, we investigate the robustness of the results of the test-retest
reproducibility analysis by running the same analysis with the Euclidean
distance between the centroid positions, defined by Eq. (14) in the main text,
to match and compare the set of centroid positions between pairs of sessions.
We also use the Euclidean distance to measure the discrepancy in terms of the
centroid position between the two optimally matched sessions. We show the $p$
values for the permutation test for the three clustering methods, five
discrepancy measures for the state-transition dynamics, and
$K\in\\{2,\ldots,10\\}$ in Tables S1, S2, and Table S3 for the DMN extracted
from the MSC data, the whole-brain network extracted from the MSC data, and
the whole-brain network extracted from the HCP data, respectively. In the
tables, $p=0$ indicates that the deviation of the ND value for the original
session-to-session comparisons from $1$ is larger than the deviation of the ND
value for all the $10^{4}$ randomized comparisons from $1$, corresponding to
$p<0.5\times 10^{-4}$. We find that the permutation test results obtained with
the Euclidean distance are similar to those obtained with the cosine distance
shown in the main text.

Quantitatively, for the DMN extracted from the MSC data, most $p$ values
(i.e., 128 out of the 135 comparisons; 94.81%) were smaller than 0.05
(indicated by *), and 111 of them (i.e., 83.07% of the 135 comparisons)
remained significant after the Bonferroni correction (indicated by ***;
equivalent to $p<0.00037$, uncorrected). Moreover, we obtained $p<0.001$,
indicated with **, in 113 out of the 135 comparisons (i.e., 83.15%). We also
found that the number of significant $p$ values (with $p<0.05$) for the
K-means (45 out of 45) and bisecting K-means (44 out of 45), which differ only
by $1$ from each other, was somewhat larger than that for the TAAHC (42 out of
45). After the Bonferroni correction, the number of the significant $p$ values
remained larger for the K-means (43 out of 45) than for bisecting K-means (36
out of 45), which is larger than for the TAAHC (32 out of 45).

Table S1: Results of the permutation test for the DMN extracted from the MSC data when we used the Euclidean distance to calculate the distance between the centroids of clusters. The entries with 0 mean that $p<0.5\cdot 10^{-4}$, uncorrected. *: $p<0.05$, uncorrected; **: $p<0.001$, uncorrected; ***: $p<0.05$, Bonferroni corrected (which is equivalent to $p<0.00037$, uncorrected). We remark that “Centroid” and “Transition prob.” abbreviate the centroid position and transition probability matrix, respectively.  |  $K$ |  Centroid |  Coverage |  Frequency |  Lifespan |  $P$
---|---|---|---|---|---|---
K-means |  2 |  $0^{***}$ |  $0.0186^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0.0014^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
TAAHC |  2 |  $0^{***}$ |  $0.5528$ |  $0.0035^{*}$ |  $0.0036^{*}$ |  $0.002^{*}$
|  3 |  $0^{***}$ |  $0.1027$ |  $0^{***}$ |  $0.0026^{*}$ |  $0.0001^{***}$
|  4 |  $0^{***}$ |  $0.0435$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0.0842$ |  $0.0001^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.0076^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  0.5548 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0.0718^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0.0434^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.2352$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
Bisecting K-means |  2 |  $0^{***}$ |  $0.0287^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0.0001^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  0.0584 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0.0419^{*}$ |  $0.0003^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.018^{*}$ |  $0.0001^{***}$ |  $0.0004^{**}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.0005^{**}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0.0037^{*}$ |  $0.0002^{***}$ |  $0.0001^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0.0072^{*}$ |  $0.0001^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.0067^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
Table S2: Results of the permutation test for the whole-brain network extracted from the MSC data when we used the Euclidean distance to calculate the distance between the centroids of clusters. The entries with 0 mean that $p<0.5\cdot 10^{-4}$, uncorrected. *: $p<0.05$, uncorrected; **: $p<0.001$, uncorrected; ***: $p<0.05$, Bonferroni corrected. |  $K$ |  Centroid |  Coverage |  Frequency |  Lifespan |  $P$
---|---|---|---|---|---|---
K-means |  2 |  $0^{***}$ |  $0.0368^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0.0402^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0.0353^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.0007^{**}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.0016^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0.0004^{**}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.0003^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
TAAHC |  2 |  $0^{***}$ |  $0.0776$ |  $0.0001^{***}$ |  $0^{***}$ |  $0.0001^{***}$
|  3 |  $0^{***}$ |  $0.0463^{*}$ |  $0^{***}$ |  $0.0003^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  0.0940 |  $0^{***}$ |  $0.0011^{*}$ |  $0.0004^{**}$
|  5 |  $0^{***}$ |  $0.0813^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.7868$ |  $0^{***}$ |  $0.0011^{*}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.8686$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  0.0875 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0.0403^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.036^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
Bisecting K-means |  2 |  $0^{***}$ |  $0.0303^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0.2482$ |  $0^{***}$ |  $0.0041^{*}$ |  $0.0174^{*}$
|  4 |  $0^{***}$ |  $0.0576$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  0.2599 |  $0.0002^{***}$ |  $0.0003^{***}$ |  $0.0483^{*}$
|  6 |  $0^{***}$ |  $0.0658$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.6427$ |  $0.0002^{***}$ |  $0^{***}$ |  $0.0001^{***}$
|  8 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  0.0905 |  $0.0018^{*}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.0031^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
Table S3: Results of the permutation test for the whole-brain network extracted from the HCP data when we used the Euclidean distance to calculate the distance between the centroids of clusters. The entries with 0 mean that $p<0.5\cdot 10^{-4}$, uncorrected. *: $p<0.05$, uncorrected; **: $p<0.001$, uncorrected; ***: $p<0.05$, Bonferroni corrected. |  $K$ |  Centroid |  Coverage |  Frequency |  Lifespan |  $P$
---|---|---|---|---|---|---
K-means |  2 |  $0^{***}$ |  $0.0023^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
TAAHC |  2 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  0.0952 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0.0159^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0.1318$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.0191^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.0096^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0.0007^{**}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0.0003^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0.0008^{**}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
Bisecting K-means |  2 |  $0^{***}$ |  $0.0016^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  3 |  $0^{***}$ |  $0.013^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  4 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  5 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  6 |  $0^{***}$ |  $0.0001^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  7 |  $0^{***}$ |  $0.0118^{*}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  8 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  9 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$
|  10 |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$ |  $0^{***}$