## Abstract

Functional connectivity (FC) describes the statistical dependence between neuronal populations or brain regions in resting-state fMRI studies and is commonly estimated as the Pearson correlation of time courses. Clustering or community detection reveals densely coupled sets of regions constituting resting-state networks or functional systems. These systems manifest most clearly when FC is sampled over longer epochs but appear to fluctuate on shorter timescales. Here, we propose a new approach to reveal temporal fluctuations in neuronal time series. Unwrapping FC signal correlations yields pairwise co-fluctuation time series, one for each node pair or edge, and allows tracking of fine-scale dynamics across the network. Co-fluctuations partition the network, at each time step, into exactly two communities. Sampled over time, the overlay of these bipartitions, a binary decomposition of the original time series, very closely approximates functional connectivity. Bipartitions exhibit characteristic spatiotemporal patterns that are reproducible across participants and imaging runs, capture individual differences, and disclose fine-scale temporal expression of functional systems. Our findings document that functional systems appear transiently and intermittently, and that FC results from the overlay of many variable instances of system expression. Potential applications of this decomposition of functional connectivity into a set of binary patterns are discussed.

## Author Summary

Numerous studies of functional connectivity have revealed densely coupled sets of brain regions corresponding to resting-state networks or functional systems. Prior work suggests that functional connectivity fluctuates over time. Here, we extend those studies by suggesting that functional connectivity can be decomposed into a set of momentary network states, with each one partitioning the network into exactly two clusters or communities. We show that these bipartitions exhibit characteristic spatiotemporal patterns that are reproducible across participants and imaging runs, and can capture individual differences. Our decomposition approach discloses fine-scale dynamics of functional systems, and reveals that functional systems coalesce and dissolve at different times and on fast timescales. Numerous applications and extensions of the approach are discussed.

## INTRODUCTION

Modern network neuroscience conceptualizes the brain as an interconnected dynamic multiscale system (Bassett & Sporns, 2017; Betzel & Bassett, 2017; Bullmore & Sporns, 2009). At the level of the whole brain, anatomical projections between brain regions shape spontaneous dynamics and constrain the brain’s momentary responses to changes in input, internal state, and environmental demand (Honey et al., 2009; Suárez, Markello, Betzel, & Misic, 2020). The resulting statistical dependencies among regional time courses are generally described as “functional connectivity,” quantified with a variety of bivariate metrics (Buckner, Krienen, & Yeo, 2013; Friston, 2011). In extant fMRI research, the Pearson correlation of blood oxygenation level dependent (BOLD) time courses remains in wide use, generally applied to long epochs of resting or task-evoked responses. The resulting correlation matrix, representing a functional network (Power et al., 2011), or “functional connectome” (Biswal et al., 2010), provides a summary representation of the system’s pairwise dependencies.

Functional connectivity, measured during the resting state, exhibits highly consistent patterns across imaging sessions (Horien, Shen, Scheinost, & Constable, 2019), participant cohorts (Dadi et al., 2019), and parcellations (Arslan et al., 2018), while also expressing individual differences (Marek et al., 2019), state-dependent changes (Betzel et al., 2020), and genetic associations (Demeter et al., 2020). Among its characteristic network features is community structure, the presence of reproducible modules consisting of regions that are internally densely coupled, reflecting their coherent and correlated activity over time. These intrinsic connectivity networks (Damoiseaux et al., 2006), or resting-state networks (RSNs), have become enshrined in the cognitive neuroscience literature, providing a fundamental taxonomy and topographic reference frame for mapping brain/behavior relations (Ito et al., 2017; Uddin, Yeo, & Spreng, 2019). Canonical sets of RSNs have been proposed (Power et al., 2011; Yeo et al., 2011), and their consistent spatial layout has been shown to reflect patterns of coactivation in task-driven fMRI activation studies (Laird et al., 2011). As internally coherent, coactivated, co-fluctuating systems they may be taken to represent building blocks of the brain’s cognitive architecture that supports specialized brain function. RSNs are not, however, sharply delineated. As has been noted in early mapping studies (Fox et al., 2005), and later revealed with data-driven community detection and clustering approaches (Power et al., 2011; Yeo et al., 2011), functional connectivity exhibits communities at multiple spatial scales, arranged in an overlapping nested hierarchy (Akiki & Abdallah, 2019; Doucet et al., 2011). Furthermore, most cognitive processes do not occur within single RSNs, and indeed may require breaking modular boundaries and dynamic reconfiguration of neural resources, including network nodes and edges (Alavash, Tune, & Obleser, 2019; Braun et al., 2015; Petersen & Sporns, 2015).

Functional systems or RSNs manifest in long-time samples of resting brain activity; indeed, their reproducibility across imaging sessions sharply increases with the length of time samples, leveling off at timescales of tens of minutes (Gordon et al., 2017). This raises the question of whether RSNs manifest only on longer timescales or whether they also “exist” at shorter timescales. Recent studies of time-varying functional connectivity (tvFC; Heitmann & Breakspear, 2018; Kucyi, Tambini, Sadaghiani, Keilholz, & Cohen, 2018; Lurie et al., 2020) have addressed the issue, approaching fine temporal structure and dynamics of FC through the use of shorter data samples, such as sliding windows or instantaneous coactivation patterns that result in temporally ordered sequences of functional networks and network states (Allen et al., 2014; Faghiri, Iraji, Damaraju, Turner, & Calhoun, 2020; Liu & Dyun, 2013; Park, Friston, Pae, Park, & Razi, 2018; Preti, Bolton, & Van De Ville, 2017; Shakil, Lee, & Keilholz, 2016). These studies have provided evidence for significant fluctuations of functional connections and network communities on timescales of tens of seconds to minutes (Grandjean et al., 2017; Hilger, Fukushima, Sporns, & Fiebach, 2020; Liao, Cao, Xia, & He, 2017; Ligeois et al., 2019; Vohryzek, Deco, Cessac, Kringelbach, & Cabral, 2020). These fluctuations involve a shifting balance between segregated (high modularity) and integrated (low modularity) states (Betzel, Fukushima, He, Zuo, & Sporns, 2016; Zalesky, Fornito, Cocchi, Gollo, & Breakspear, 2014), with episodes of high modularity exhibiting consistent topology across time (Fukushima et al., 2018) and subject to modulation by internal state, task performance, or behavior (Cohen, 2018; Shine et al., 2016).

Recently, we suggested a new approach to functional connectivity, by focusing on the dynamics and networks formed by “edge time series” (Esfahlani et al., 2020; Faskowitz, Esfahlani, Jo, Sporns, & Betzel, 2020; Jo, Faskowitz, Esfahlani, Sporns, & Betzel, 2020). The approach unwraps time-averaged FC into time series of co-fluctuating signals on network edges resolved at the timescale of single frames in MRI acquisition, thus allowing inspection of network dynamics at fine timescales. Here we build on this approach and show that a simple proxy for the resulting framewise community structure, expressed as a set of bipartitions of the network into two positively co-fluctuating ensembles of nodes, represents a compact decomposition of the functional connectivity. Examining the patterns and frequencies of these bipartitions allows addressing several issues related to FC dynamics. How well do bipartitions, sampled at fine-scale temporal resolution, represent “classic” system-level architecture as derived from long-time FC? How does the community structure of single frames combine into the community structure of FC? Do systems, as coherent blocks, manifest continuously or do they appear transiently at short timescales? Do systems differ in their patterns of “functional expression” across time?

## RESULTS

### Extraction of Bipartitions From Time Series

This expository section introduces the basic constructs employed in this study (Figure 1); for more detail see the Methods section.

Starting from node time series (BOLD activations), the cross-correlation between
each pair of nodes defines their linear statistical dependence (Figure 1A). The correlations of all node
pairs within a given system are that system’s functional connectivity.
Employing a standard definition of the cross-correlation, the average of the
products of the standard scores of the two variables, yields scalar correlation
estimates. Omitting the averaging step retains the summands, corresponding to a
temporal unwrapping of the scalar correlation estimates into vectors (time
series) along each edge (corresponding to the link between a node pair). These
edge time series represent co-fluctuations of node pairs, which are positive when the sign of
the two nodes’ signal amplitude agrees, and negative otherwise. The
average of these edge time series is equivalent to the value of the
corresponding correlation (functional connectivity) and, when computed across
all edges, is equivalent to the FC matrix (Figure 1B). A useful summary metric aggregates the amplitudes of all edge
co-fluctuations, computed as the square root of the sum of their squared values
(root sum square), here denoted *RSS*. High *RSS* values indicate that node signals strongly agree/disagree at a given point in
time.

Removing amplitudes and retaining only the sign of co-fluctuation along edges
naturally partitions the network into exactly two sets of nodes (Figure 1C), one set comprising nodes with
positive z-scores and a complementary set comprising the remaining nodes with
negative z-scores. This is equivalent to thresholding each frame’s node
vector at *z* = 0. These two sets of nodes internally
co-fluctuate positively and exhibit negative co-fluctuations between them, thus
defining a bipartition of the network. Reversing the sign of BOLD amplitudes
will retain the exact same co-fluctuation pattern and bipartition; we will
therefore disregard the signs of z-scored node amplitudes in further analysis.
Note also that applying the *z* = 0 threshold, while inherent to
the computation of FC from edge time series, should not imply functional
activation of nodes above *z* = 0; it merely indicates that
regions are active above or below their own mean.

Bipartitions divide the network into exactly two communities, and, over all time frames, these community assignments can be combined into a co-assignment or agreement matrix. In network science, agreement matrices are often used to represent graded assessments of community affiliation (also called co-classification or co-assignment), for example in consensus clustering (Lancichinetti & Fortunato, 2012) and multiresolution community detection (Jeub, Sporns, & Fortunato, 2018). In general, the agreement matrix expresses the frequency with which each node pair is grouped into the same community across many partitions. Here, we calculate the agreement of many bipartitions across many time points.

Bipartitions, as special cases of partitions that bisect the network into two
communities, are described by a binary node vector of community assignments. The
similarity between two such vectors can be measured with several distance
metrics such as the Jaccard distance, the cosine similarity, the variation of
information, or the mutual information. Here, we adopt mutual information
(*MI*) as the principal metric used for assessing similarity
between bipartitions. Other metrics are highly correlated with *MI*, and their application gives qualitatively similar
results to those reported in this article.

Variations of the bipartition approach are possible, and a few are explored as
part of this study. For example, the zero-threshold dividing each frame into two
sets of nodes based on the sign of their z-scored time courses may be modified
by adopting an arbitrary threshold *θ*. Another approach
is to define two thresholds +*θ* and
−*θ* that separate highly positive and highly
negative activations from activations near the temporal mean, thus yielding
tripartitions.

### Bipartitions Are Strongly Related to Functional Connectivity

All analyses reported in this article have been carried out on resting-state fMRI scans acquired in a cohort of 95 participants, a quality-controlled subset of the “100 unrelated” Human Connectome Project (Glasser et al., 2013) cohort. After preprocessing and nuisance regression, each one of the four separate runs was composed of 1,100 frames (TR = 720 ms, total length 792 s). BOLD time courses from cerebral cortex were parcellated into 200 nodes according to a standard template (Schaefer et al., 2018). Some variations of MRI preprocessing have been explored and are discussed below, including a second parcellation scheme into a finer set of 300 nodes and an alternative nuisance regression strategy that omits global signal regression (referred to as “non-GSR data”). For details on participants, scanning, and fMRI processing see the Methods section. To allow division of brain regions into a set of functional systems, each network node was assigned to one of seven canonical RSNs (Yeo et al., 2011), comprising the visual (VIS), somatomotor (SOM), dorsal attention (DAN), ventral attention (VAN), limbic (LIM), frontoparietal (FP), and default mode (DMN) systems (Figure S1 in the Supporting Information).

Classic FC is equal to the mean over all frames (time points) of the edge time
series (Figure 2A). Edge time series are
converted to binary form by applying a threshold based on the sign of the
momentary co-fluctuation, an operation that results in a series of bipartitions
(Figure 2B). The agreement matrix
constructed from these bipartitions is highly correlated with the corresponding
FC matrix ($r^$ = 0.964 ± 0.008, 95 participants, one run). The strong correlation
between this bipartition overlay and traditional that FC is robust against
different choices of node parcellation and fMRI preprocessing. Figure S2 shows consistently strong
similarity between FC and agreement matrix for a finer nodal parcellation (300
nodes) and for time series data omitting global signal regression.
Interestingly, for both variants of preprocessing the agreement matrix contains
negative entries, representing node pairs whose co-assignment into the same
module was below the level predicted by chance. In global-signal-regressed data,
these entries strongly overlap with negative functional connectivity. Variants
of the bipartition approach also yield high matches of agreement and FC
matrices. Adopting an arbitrary (nonzero) threshold *θ* to
create bipartitions, or adopting an approach using tripartitions, results in
close approximation of FC over a wide range of the *θ* parameter (Figure S3).

The set of bipartitions is a binary decomposition of functional connectivity. The
characteristic patterning of FC is constructed from the specific spatiotemporal
patterns of its constituent bipartitions, as shown in Figure 3. Subsampling randomly chosen bipartitions gradually
approximates FC, with even modest proportions of frames (around 10%)
resulting in a very close match with the full-length FC estimate (Figure S4A). When varying run length and
using all frames, the quality of the match between agreement and FC matrices
remains high even when runs are short (Figure
S4B). Prior work (Esfahlani et al., 2020) noted that selecting edge time series frames based on their
rankings in *RSS* magnitude approximates FC more quickly when
frames are ranked from high to low *RSS* amplitudes, as opposed
to ranking them from low to high. Bipartitions behave very similarly (Figure S4C). This effect persists when
accounting for the autocorrelation structure (temporal adjacency) of the
selected frames (Figure S4D). The level
to which bipartitions approximate FC is unrelated to framewise head motion.
“Scrubbing” (removing) high-motion frames (retaining only the
frames below the 90th percentile of the framewise displacement, per participant,
per run) does not significantly affect the match between agreement and FC
matrices (Figure S5).

Representing the fMRI time series as a series of bipartitions allows computing
their pairwise similarity (quantified as mutual information) across time. Figure 3A displays an example of such a
similarity matrix for a single participant and a single run. Notably, some
instances of bipartitions recur throughout the run, as indicated by strongly
positive MI between remote time points (off-diagonal entries in the matrix
plot). Reordering frames by *RSS* magnitude on each run, followed
by averaging over all participants, reveals that high similarity of bipartitions
is largely restricted to episodes when *RSS* amplitudes are near
maximal (Figure 3B). The bipartition
similarity between adjacent time points (pairwise MI) is correlated with *RSS* amplitude (Figure 3C shows data from a representative participant; *ρ* = 0.502, *p* =
10^{−71}; $\rho ^$ = 0.494 ± 0.049, 95 participants, one run). Lower values of pairwise MI
occur when *RSS* amplitudes are small, and higher pairwise MI
occurs predominantly when *RSS* amplitudes are large. This
relationship indicates that the community structure expressed in framewise
bipartitions is more stable (changes less) when overall co-fluctuations, across
the entire network, are large. These time points correspond to moments when BOLD
time series (and hence co-fluctuations), on average, exhibit larger amplitudes,
that is, are farther from their zero mean. Different nodes switch at different
rates (Figure 3D), with several DMN regions
(parcels DefA_IPL_1 and DefA_PCC_1, both hemispheres) and VAN regions (parcel
SalVentA_ParOper_1, right hemisphere) remaining most stably associated with
their host communities during high-amplitude epochs.

Principal component analysis (PCA), applied to the set of bipartitions extracted
from each participant’s BOLD time course, yields a small number of
principal components (PCs) that account for significant portions of the observed
variance and exhibit consistent topography across participants (Figure 4A). The largest PC (PC1), on average,
accounts for approximately 12% of the variance (12.772 ± 2.272,
range 20.393 to 8.970, 95 participants, one run). Averaged across participants
and projected onto the cortical surface, the PC1 pattern corresponds to a mode
that splits the brain into two co-fluctuating ensembles comprising most regions
belonging to the VIS, SOM, DAN, and VAN systems on one side versus most regions
belonging to the LIM, FP, and DMN systems on the other (Figure 4B). The temporal expression of PC1 is positively
correlated with *RSS* amplitude (Figure 4C), indicating that the connectivity mode inscribed in PC1
is most strongly expressed at time points with high-amplitude network-wide
co-fluctuations. This correlation with *RSS* is obtained despite
discarding amplitude information when creating binary partitions. The PC1 as
derived from sets of bipartitions is related to several other more familiar
constructs (Figure 4D). It is equivalent to
the principal eigenvector of FC (or, more precisely, its corresponding
covariance matrix), and thus also exhibits strong resemblance to connectivity
“gradients” (Margulies et al., 2016). We derived principal components of the affinity matrix
derived from FC (equivalent to the principal FC eigenmode) as in the example
shown in Figure 4D. The resulting pattern
is very highly correlated with the PC1 derived from bipartitions. Furthermore,
the bipartition PC1 pattern closely resembles the first principal component of
bipartitions identified by applying the Louvain modularity maximization
algorithm to the long-time-averaged FC matrix (Figure 4D). These strong relationships indicate that the set of
bipartitions encapsulates characteristic features of the FC matrix, including
its eigenmodes and community structure.

Collectively, these results demonstrate a strong relationship between bipartitions and traditional functional connectivity as expressed in the FC matrix. The set of bipartitions derived from the BOLD time series reflects several important spatial and topographic features of FC, while also disclosing its fine temporal structure.

### Bipartitions Map Onto Basis Sets of Templates

Bipartitions divide the network, at each point in time, into exactly two communities. These two communities are often approximately equal in size, with only 5% comprising node sets that have fewer than 70 (out of 200) members. This fact begs the question of how these large communities relate to canonical subdivisions of the brain into much more compact functional systems, the largest of which (the DMN in the 200-node parcellation) comprising 46 nodes. One way to address this question is to compare each of the empirically observed bipartitions with a standard or basis set of templates that split the brain into bipartitions defined along the boundaries of canonical functional systems (Figure 5A). The basis set used here comprises 7 templates that divide 7 canonical RSNs (Yeo et al., 2011) into 1 versus 6 networks, 21 templates that divide them into 2 versus 5 networks, and 35 templates that divide them into 3 versus 4 networks, for a total of 63 such templates. Since these templates are drawn along RSN boundaries (which themselves are defined based on their coherent co-fluctuations over long timescales), one would expect that bipartitions observed at each frame will at least partially align with the boundaries of the 7 systems as captured in the 63-template basis set.

Comparison of templates with observed bipartitions over time allows tracking of
several metrics: (a) the similarity (mutual information) of each bipartition
with each basis set template; (b) the identification of the single basis set
template that most closely resembles the observed bipartition; and (c) computing
which of these best-matched templates occur most frequently and which correlate
most strongly with framewise measures such as *RSS* amplitude. Figure 5B shows examples of three MI
time courses for three examples of templates (cf. Figure 5A), one each that divides the network into 1 + 6, 2 + 5, and
3 + 4 systems. The full set of 63 MI time courses represent how well each
observed bipartition resembles each of the 63 basis set templates and may be
interpreted as an index of how strongly a given template is realized at a given
point in time. Selecting, at each time frame, the template for which the MI is
maximal allows representing the sequence of highly variable bipartitions as a
sequence of integers, each representing the single best match (highest MI) out
of the 63 templates. Figure S6 provides
examples of observed bipartitions and their best matches in the template set
determined by maximal MI, for three example templates.

For each participant and imaging run, templates can be ordered by their median
frequency, based on the number of times they were selected as the best match for
the observed bipartitions (Figure 5C). Once
a single best-matching basis set template is assigned to each frame, their
occurrence can be compared against *RSS* amplitude (Figure 5D). The most frequently observed
basis set template (template 63) most strongly correlates with framewise *RSS*, indicating that it is predominantly expressed when
BOLD signals and their co-fluctuation patterns exhibit high amplitudes. Note
that the template 63 pattern strongly resembles the PC1 extracted from observed
bipartitions (cf. Figure 5A). Qualitatively
similar rankings of basis set templates and correlations with *RSS* amplitude are obtained for a finer node parcellation
and for non-GSR data (Figure S7).

Finally, the best-matching template set represents a highly compressed set of features of the framewise decomposition, specific to each imaging session and to each participant. Discarding the temporal ordering of the templates, which is immaterial for computing or reconstructing FC, results in a string of 63 numbers encoding a frequency spectrum. The agreement matrix of the template set, as encoded in the 63-element vector, closely matches the down-sampled system-wise FC (Figure S7). This suggests that the frequency of the template basis set may represent a highly compact encoding of the empirically observed FC.

### Bipartitions Contain Signatures of Individual Differences

If templates represent an efficient encoding of the framewise decomposition of the full-length FC time course, can features of this encoding be useful for identifying individual differences in FC? One way to approach this question is to compare individual variations in template frequencies across multiple resting-state runs (4 runs yielding 6 possible pairwise comparisons: runs 1 vs. 2, 1 vs. 3, 1 vs. 4, 2 vs. 3, 2 vs. 4, 3 vs. 4). Results show that frequencies of a subset of templates are significantly correlated when comparing pairs of runs within participants (Figure 6A). This suggests that, at least for some template classes, the level of expression during resting state is moderately conserved across individuals.

A second way to address the question is to pursue connectional
“fingerprinting” (Finn et al., 2015) by comparing patterns of functional connectivity across
individuals and across runs. Similarity of FC across the cohort is assessed by
computing the differential identifiability and accuracy (see the Methods section), both computed for all 6 pairwise
combinations of runs and then averaged. Using all 1,100 frames per run to
compute FC yields *I*_{diff} = 12.97 and
accuracy = 68.58 (mean of 6 comparisons). Substituting the corresponding
agreement matrix yields *I*_{diff} =
12.32 and accuracy = 65.75. Framewise decomposition of FC time series offers the
opportunity to select subsets of frames, derive the means of their edge time
series and the agreement of their bipartitions, and then compute *I*_{diff} and accuracy for these
samples.

Frames may be selected by applying a wide range of criteria. Here, we selected a
criterion that was designed to return frames with bipartitions that were closest
to each of the 63 templates in the basis set (Figure 5A). The criterion was computed as follows. For each
participant and run, we previously calculated time series of MI between the
observed bipartition and each of 63 templates (see Figure 5B). At each time step, the set of 63 MI values was converted
to z-scores. For each participant, each run, and all templates we then selected
the top 10% of frames with the highest z-scores, comprising a subset of
frames with bipartitions that were highly similar to a given target template.
Identifiability and accuracy computed from the selected subset are compared
against multiple (20) instantiations of a null model where the selected time
points (for each participant, each run) were shifted by a random value. Figures 6B and 6C show *I*_{diff} and
accuracy for each of the 63 target templates. Results suggest that several
target templates significantly outperform their respective nulls. Highest levels
of identifiability and accuracy are attained when subsets of frames are selected
based on their similarity to template 4, representing high co-fluctuation within
the VAN system, and template 24, representing high co-fluctuations within the
VAN and FP systems.

Figure 6D shows normalized (z-scored)
similarity matrices, displayed as the mean of 6 matrices of Pearson correlations
of subsampled mean edge time series computed across all pairs of participants.
Values on the main diagonal refer to self-similarities of participants across
runs, while off-diagonal values refer to similarities between different
participants across scans. Normalization suggests that increases in *I*_{diff} and accuracy for target
template 4 relative to target template 63 are due to increased self-similarity
(cf. Finn et al., 2015). Figure 6E shows the between-participant
correlations of the mean agreement matrix (means over all 4 runs) for selected
subsets of frames. Selecting subsets of frames based on target template 4
results in lower between-participant correlations ($r^$ = 0.562) than selecting them based on template 63 ($r^$ = 0.835; corresponding means for edge time series patterns: $r^$ = 0.601 and $r^$ = 0.828; mean correlations when using all 1,100 frames: $r^$ = 0.740 for agreement and $r^$ = 0.752 for edge time series). Figure 6F displays the average agreement matrix (across all participants and runs) for
target templates 4 and 63.

Collectively, these findings suggest that framewise decomposition and template classification may exhibit some level of stable individual differences and may be useful for enhancing participant-specific network patterns. Regarding the latter point, it appears that selecting moments in time when VAN and/or FP are strongly co-fluctuating results in network patterns that express enhanced levels of individual differences.

### Bipartitions Are Constrained by Functional Connectivity

So far, findings indicate that the set of bipartitions observed during single resting-state fMRI runs closely approximates FC (Figure 2) and exhibits characteristic spatiotemporal patterns (Figures 3, 4, 5, and 6). Working backwards from a given FC matrix, we can ask to what extent the long-term pattern constrains the set of underlying fine-scale bipartitions from which it is composed. Obviously, many different sets of bipartitions (many different sets of time courses) can yield identical FC. To what extent are sets of bipartitions free to vary once their final superposition in FC is fixed?

An optimization approach, searching the space of all possible bipartitions, can help address this question (Figure 7). The approach adopts a variant of the Metropolis algorithm (Metropolis, Rosenbluth, Rosenbluth, Teller, & Teller, 1953) by maximizing an objective function, defined as the similarity between an empirically observed agreement matrix (which, as established above, very closely resembles FC) and an agreement matrix derived from a set of bipartitions that are subject to incremental optimization. The initial state consists of a completely random set of bipartitions that give rise to a flat agreement matrix. Then, at each subsequent iteration, a single node’s community affiliation on a single time frame (both chosen uniformly and randomly) is swapped. The objective function is recomputed after each swap, and the swap is retained if similarity is increased, subject to a simulated annealing paradigm (Kirkpatrick, Gelatt, & Vecchi, 1983) applied to ensure that the end state corresponds, as closely as possible, to a global optimum. Three different objective functions are employed, the Pearson correlation, the cosine distance, and the root-mean-square distance (additional data shown in Figures S9A and S9B), with near-identical outcomes. Applying the algorithm to data from single participants and single imaging runs succeeds in identifying sets of bipartitions that closely approximate the agreement matrix derived from the empirical BOLD time series (Figure 7A). Importantly, the optimized set of bipartitions resembles the set of observed bipartitions, as determined by comparing their respective best-matching basis set templates (Figure 7B; Figure S9C). Optimization also yields closely matching sets of bipartitions when the optimized set of bipartitions is significantly smaller than the length of the original time series (Figure S9D). For example, if the optimized set is limited to one-tenth of the length of the original time series (110 frames), optimization still converges and resulting bipartitions continue to resemble those in the observed set.

These findings suggest that the set of bipartitions encountered in the decomposition of fMRI data is highly constrained by the long-time average functional connectivity. Recall that each bipartition represents a snapshot of how momentary co-fluctuations distribute across the network, and that the total set of these snapshots exhibits significant fluctuations across time. The optimization approach suggests that these fluctuations are necessary for reconstructing long-time averages in FC, as optimized bipartitions strongly resemble the highly variable observed set.

### Expression of Canonical Systems Varies Across Time

The findings presented so far suggest that bipartitions offer an opportunity to compress time courses into discrete feature sets that retain long-time characteristics of FC while also disclosing fine-scale dynamics. A complementary approach to extract fine-scale network states is possible, as explored in this final section. The expression of individual functional systems across time can be tracked directly, by examining co-fluctuation patterns at fine-scale temporal resolution. The mean co-fluctuation of functional systems can be computed across all 7 × 7 subblocks (each system and each system interaction), yielding 28 unique time series. An example is shown in Figure 8A. The temporal averages of these time series are identical to the corresponding down-sampled 7 × 7 functional connectivity matrix (cf. Figure S8). On each time step, mean co-fluctuations are compared with a null distribution derived by randomly shuffling system labels and recomputing co-fluctuations (100 independent shuffles per time step). This comparison yields z-scores for each system and pairwise system interaction, where the z-score expresses how much the signal deviates from the label-reshuffling null. Discretizing these time courses by applying a z-score threshold (Figure 8B) yields discrete “network states,” with systems and between-system interactions either exceeding or failing to exceed the threshold of expression. Visual inspection of the sample time course suggests that each of the seven RSNs is significantly expressed, as indicated by exceeding the co-fluctuation z-score threshold, only intermittently, on a fraction of time points. Recall that co-fluctuations should not be taken as “mean activation time courses” as they take on positive values when participating nodes are either jointly above or jointly below their long-time z = 0 means.

Considering the above-threshold expression of each of the seven RSNs (leaving
aside their mutual interactions, and noting that strongly negative z-scores do
not occur) yields, for each point in time, a binary seven-element vector (a
total of 128 such states are possible, with between 0 and 7 RSNs expressed at a
given time). Aggregating these states (95 participants, 4 runs, 418,000 frames)
provides summary statistics on their frequency (Figure 8C). The most frequent state (occurring in approximately
13% of all frames) is one where no RSN is strongly expressed. Individual
participants range between 7.4% and 22.4%, and expression levels
are correlated across imaging runs (95 participants, $\rho ^$ = 0.375, p < 10^{−4} for five out of six pairwise
correlations; Figure S10A). The next most
frequent states predominantly include those where single RSNs are significantly
expressed, while states that involve simultaneous co-expression of multiple
systems are less frequent. Frequency ranks of states remain stable when changing
z-thresholds (Figure S10B).

As discussed above, bipartitions decompose FC into sequences of two-community assignment vectors that can be matched to templates from a basis set. As defined in this section, network states also represent sequences of discrete patterns directly derived from significant excursions of edge time series. How do these two representations relate to each other? Network states derived from system-wise expression levels partially reflect the community structure of bipartitions. Many of the states expressing one or several canonical functional systems have clear counterparts within the bipartition template set, that is, the two representations coincide in time (Figure 8D). Aggregating the bipartitions observed on each time point corresponding to the 10 most frequent network states confirms that most states map onto consistent patterns of co-fluctuation as indexed by the bipartition approach (Figure 8E). This comparison establishes a relationship between network states as defined here, through framewise averaging of co-fluctuations, and the community structure of bipartitions as defined in previous sections. Both represent compact descriptions of the dynamic expression of functional systems on fine timescales.

## DISCUSSION

Fine-scale analysis of BOLD signal co-fluctuations (edge time series) demonstrates that canonical functional systems are not expressed uniformly or stably across time. Instead, their levels of expression fluctuate significantly, as individual functional systems coalesce and dissolve, singly or in varying combinations. While found reliably and reproducibly in long-timescale FC, this appearance is the result of the overlap of many transient and fleeting manifestations. The proposed decomposition of FC into bipartitions and network states allows tracking these dynamics at fine timescales, only limited by the acquisition rate of single MRI frames. The approach complements traditional network analysis of FC, estimated on long timescales, and of time-varying FC, estimated on shorter windows or epochs.

We propose that FC can be decomposed into sets of bipartitions that map onto discrete network states. These bipartitions exhibit characteristic spatiotemporal patterns, with systems and combinations of systems expressed at different times, and in varying combinations. The most common patterns are those where none of the systems are expressed, or where systems are expressed singly and in isolation. The statistics of bipartitions and network states are reproducible across imaging runs and participants, and do not appear to depend critically on choices made in fMRI preprocessing (e.g., parcellations and global signal regression). Their patterning reflects the complex multiscale community structure of long-time FC, which has been, to this point, the primary target of functional network analysis.

Our work builds on and extends previous investigations of time-varying functional
connectivity that has provided evidence for time-dependent fluctuations in
functional connections (Chang &Glover, 2010) and network patterns and states (Allen et al., 2014; Lurie et al., 2020; Zalesky et al., 2014).
Consistent with prior studies of tvFC, our approach reveals spatiotemporal patterns
of network-wide co-fluctuations. Notably, we detect a dominant (segregated or
modular) connectivity mode that covaries with overall signal amplitudes
(*RSS*), appears intermittently over time, and exhibits
consistent topography (Betzel et al., 2016;
Fukushima et al., 2018; Shine et al., 2016). Unlike most tvFC studies, our approach
does not require defining sliding windows and hence allows tracking system dynamics
at higher temporal resolution. The decomposition of the edge time series into
discrete sets of bipartitions and/or network states offers not only a highly
compressed encoding of system dynamics but also potential new targets for analysis
and modeling of both resting and task-evoked fMRI time series data. Our approach is
related to other methods, including coactivation pattern (CAP) analysis
(Karahanoğlu & Van De Ville, 2015; Liu & Dyun, 2013;
Liu, Zhang, Chang, & Duyn, 2018) as
well as approaches for detecting transient network states that employ hidden Markov
models (Vidaurre et al., 2018; Vidaurre,
Smith, & Woolrich, 2017), dynamic
effective connectivity (Park et al., 2018),
instantaneous phase synchrony (Pedersen, Omidvarnia, Walz, Zalesky, &
Jackson, 2017; Pedersen, Omidvarnia,
Zalesky, & Jackson, 2018), patterns
of instantaneous phase-locked BOLD responses (Vohryzek et al., 2020), or filter banks (Faghiri et al., 2020). In contrast to most of these approaches, our
method is distinct in that it does not require the user to specify a seed region or
a threshold for an “event” (Allan et al., 2015; Petridou, Gaudes, Dryden, Francis, &
Gowland, 2013; Tagliazucchi, Balenzuela,
Fraiman, & Chialvo, 2012), and does
not require “windowing” or specification of models or parameters for
the inference of states and state transitions. Importantly, the temporal average of
edge time series generated by our approach is exactly equal to time-averaged FC,
making it possible to retrieve the precise contributions of individual frames to the
static correlation pattern.

The topography of the dominant principal component of the bipartitions bears strong
resemblance to the principal mode of BOLD dynamics observed during high *RSS* amplitude “events” (Esfahlani et al., 2020), as well as patterns characterized by
strong excursions (Betzel et al., 2016) or
high modularity (Fukushima et al., 2018) in
time-varying functional connectivity. Similar patterns representing a decoupling of
mainly task-positive from task-negative regions have been described and interpreted
in previous studies as a major intrinsic/extrinsic dichotomy in functional
architecture (Doucet et al., 2011; Fox et
al., 2005; Golland, Golland, Bentin,
& Malach, 2008; Zhang et al., 2019). The pattern reported here is also very
highly correlated with cortical gradients (Margulies et al., 2016), specifically those derived from
eigen-decompositions of the functional connectivity matrix. Indeed, this strong
resemblance is due to a mathematical relationship between sets of framewise
bipartitions described here (a compression of the original time series) and the
spatial patterns of FC eigenmodes. Going beyond static patterns such as gradients
(see also Faghiri, Stephen, Wang, Wilson, & Calhoun, 2019), our approach links these connectivity eigenmodes
to fluctuating levels of expression of specific functional systems at fine-scale
temporal resolution. Their relation to cognitive processes, such as ongoing thought
(Mckeown et al., 2020), is a topic for
future study.

Our findings suggest that the decomposition of FC into edge time series and bipartitions may offer new approaches for extracting network markers of individual differences. Prior work has established that long-time averaged FC contains connectional features that allow connectotyping (Miranda-Dominguez et al., 2014) or fingerprinting (Finn et al., 2015) of individuals. Subsequent studies have shown that signatures of individual variability involve specific functional systems (Finn et al., 2015), can be extracted from small sets of nodes and connections (Byrge & Kennedy, 2019), and can be enhanced by extracting subsets of principal components (Amico & Goñi, 2018). Analysis of tvFC patterns suggested that cortical fingerprints are expressed at different levels across time, with FC snapshots that are most unique across individuals carrying the most information regarding individual differences (Peña-Gómez, Avena-Koenigsberger, Sepulcre, & Sporns, 2018). Here, we build on this body of work and suggest that FC and agreement matrices derived from subsets of frames may carry enhanced signatures of individual differences. Our findings suggest that frames characterized by strong co-fluctuations within the ventral attention and frontoparietal networks carry information that is specific to participants. These systems, or their close homologs, have also been implied in previous work on FC fingerprinting (Finn et al., 2015; Peña-Gómez et al., 2018). An intriguing avenue for further study, suggested by findings on participant-specific and heritable state frequencies and transitions (Vidaurre et al., 2017), is whether decomposition of FC and selection of frames can uncover new connectivity traits related to cognition and behavior.

We employed an optimization approach to explore whether a given FC matrix can be decomposed into sets of bipartitions that differ radically from the ones that are empirically observed. Our findings suggest this is not the case. Given an observed pattern of FC, the set of framewise patterns from which it is composed, or into which it can be decomposed, is not free to vary. Instead, the statistics of these patterns appear strongly constrained by the correlation structure inscribed in long-time FC. Optimizations invariably retrieve sets of patterns that resemble those observed empirically, even though no dynamic generative model is employed. This makes it harder to dismiss the observed framewise patterns as artifactual or as massively corrupted by noise or uncertainty. Instead it appears that the fluctuating patterns of observed bipartitions are necessary, in the sense that it is difficult if not impossible to construct FC from a radically different (“stationary,” temporally smoother, less variable) set of frames. Our findings suggest the hypothesis that framewise bipartitions reflect a decomposition of the community structure of the FC matrix, unwrapped in time.

Our work opens new avenues for future research. The decomposition of BOLD time series into a set of bipartitions and/or network states represents a compression or encoding of the system’s dynamics into a much more compact feature set. Such feature sets may provide novel opportunities for mapping individual differences, relations to demographic or behavioral measures, task-rest reconfigurations, or relation to the underlying anatomy. They may also serve as tools for further analysis using machine learning or multivariate statistical mapping, including those probing the relation of brain to behavior. In addition, we note that the proposed scheme may also apply to brain data obtained with other acquisition methods, including more highly resolved recordings of neuronal populations or individual neurons. The decomposition of FC into framewise contributions allows us to selectively recombine subsets of frames to get different patterns of FC. It might be possible to select specific subsets of frames/templates to amplify a brain/behavior correlation. Another possibility for future research lies in exploring alternative basis sets of templates. Here, many of our findings related to bipartitions are expressed by reference to a specific choice involving a 63-template basis set. This choice was motivated by a specific research agenda related to the expression of canonical functional systems for which a nodal partition, and hence basis set, are defined. However, many alternative choices of templates are possible, depending on specific working hypotheses. Templates could also be defined in a more data-driven way, as centroids derived by clustering techniques, or from modular partitions of the empirically measured FC matrix. Finally, alternative decomposition strategies and edge metrics may be pursued. Substituting FC by the covariance matrix, a framewise decomposition could retain information on instantaneous BOLD signal variance, a measure that has proven valuable in prior work (Garrett, Kovacevic, McIntosh, & Grady, 2010). Substituting time-lagged correlation, partial correlation, or covariance, or information-theoretic measures such as mutual information, may provide additional opportunities for decomposing long-time averages for node pair interactions into framewise contributions that disclose temporal patterns and fluctuations.

Limitations of the approach should be noted. As is the case with all studies employing functional neuroimaging, the present work inherits most drawbacks of fMRI methodology including its limited temporal and spatial resolution, the indirect link to underlying neural activity, and measurement noise and statistical biases. Following good practices in data preprocessing, the use of multiple data sources, and cautious interpretation of findings can at least partially guard against these limitations. It should also be noted that the basic methodological framework transcends the limitations of fMRI as its mathematical and algorithmic core applies to all time-dependent data sources regardless of origin. Future work should aim to reproduce and refine the proposed feature sets, spatiotemporal patterns, and statistics on system expression by leveraging new data sources and participant cohorts. Interventional studies and multimodal experimentation are needed to identify putative neurobiological mechanisms that underpin or drive temporal fluctuations.

Other limitations relate more directly to the analytic approach. The calculation of cross-correlations as measures of BOLD functional connectivity implies z-scoring of the original time series, which effectively normalizes and centers all time series. If means are nonstationary, for example due to task-induced regional activations, more sophisticated methods that correct for such transient responses are needed (e.g., Cole et al., 2019). Finally, as has been pointed out in the case of “windowing” approaches in tvFC, finer temporal scales (shorter time segments) may increase the noisiness of each estimate of FC (Leonardi & Van De Ville, 2015; Zalesky & Breakspear, 2015). In our approach, no windowing is applied, and each co-fluctuation pattern is directly derived from framewise BOLD responses. It is worth pointing out that long-time FC and noisy BOLD patterns are intertwined. Adding noise to framewise BOLD signals would imply adding noise to the FC as well, and vice versa. Because the decomposition into frames and bipartitions is mathematically exact, it is not subject to noisy estimation.

In conclusion, fine-scale temporal fluctuations in the community structure of resting brain activity suggest that the brain’s functional systems express only transiently, intermittently, and infrequently across time. Their robust manifestation over long timescales results from the superposition of large numbers of spatially distinct and temporally variable patterns. Novel insights and applications may result from the proposed decomposition of brain dynamics into network bipartitions and states.

## MATERIALS AND METHODS

### Dataset and fMRI Preprocessing

All analyses reported in this article were carried out on data originally collected by the Human Connectome Project (HCP; Van Essen et al., 2013), specifically resting-state fMRI data from 100 unrelated adult participants (54% female, mean age = 29.11 +/−3.67, age range = 22–36). The study was approved by the Washington University Institutional Review Board and informed consent was obtained from all subjects. Participants underwent four 15-min resting-state fMRI scans (here referred to as four imaging runs) spread out over a 2-day span. For a full description of the imaging parameters and image preprocessing, see Glasser et al. (2013). Briefly, data were acquired with a gradient-echo EPI sequence (run duration = 14:33 min, TR = 720 ms, TE = 33.1 ms, flip angle = 52, 2-mm isotropic voxel resolution, multiband factor = 8). Participants were instructed to keep their eyes open and fixate on a cross. Images were collected on a 3T Siemens Connectome Skyra with a 32-channel head coil. Subjects were considered for data exclusion based on the mean and mean absolute deviation of the relative root-mean-square motion across either four resting-state MRI scans (file: Movement_RelativeRMS.txt) or one diffusion MRI scan (file: eddy_unwarped_images.eddy_movement_rms), resulting in four summary motion measures. If a subject exceeded 1.5 times the interquartile range (in the adverse direction) of the measurement distribution in two or more of these measures, the subject was excluded. These exclusion criteria were established before the current study. Four subjects were excluded based on these criteria. One subject was excluded for software error during diffusion MRI processing. Even though diffusion MRI was not part of the present study, this subset was created to include subjects with adequate resting-state and diffusion data for future analysis. The remaining subset of 95 subjects have the following demographic characteristics: 56% female, mean age = 29.29 +/−3.66, age range = 22–36. Finally, we note here that we defined framewise displacement as the relative root-mean-square motion (file: Movement_RelativeRMS.txt), which was computed with the FSL function rmsdiff via the HCP pipelines. We used this information to censor the resting-state scans at a frame-by-frame level in a supplementary analysis.

HCP data were minimally preprocessed as described in Glasser et al. (2013). Briefly, data were corrected for gradient distortion, susceptibility distortion, and motion, and then aligned to a corresponding T1-weighted (T1w) image with one spline interpolation step. This volume was further corrected for intensity bias and normalized to a mean of 10,000. This volume was then projected onto the 32k fs LR mesh, excluding outliers, and aligned to a common space using a multimodal surface registration (Robinson et al., 2014).

A functional parcellation designed to optimize both local gradient and global similarity measures of the fMRI signal (Schaefer et al., 2018; Schaefer200) was used to define 200 regions (parcels or nodes) of the cerebral cortex. These nodes can be mapped to a set of canonical functional networks (Schaefer et al., 2018; Yeo et al., 2011); in the current study we adopt a mapping to seven canonical networks that comprise the visual (VIS), somatomotor (SOM), dorsal attention (DAN), ventral attention (VAN), limbic (LIM), frontoparietal (FP), and default mode (DMN) systems (cf. Figure S1). For HCP data, the Schaefer200 is openly available in 32k fs LR space as a cifti file. A second processing variant used a finer parcellation into 300 nodes (Schaefer300) following the same basic procedure.

We employed two variants in preprocessing to explore the robustness of main findings reported in this article. All analyses were first carried out on data processed with the inclusion of global signal regression. For this strategy, the mean BOLD signal for each cortical node was linearly detrended, band-pass filtered (0.008–0.08 Hz; Parkes, Fulcher, Ycel, & Fornito, 2018), confound regressed, and standardized using Nilearn’s signal.clean, which removes confounds orthogonally to the temporal filters (Lindquist, Geuter, Wager, & Cao, 2019). The confound regression employed (Satterthwaite et al., 2013) included six motion estimates, time series of the mean CSF, mean WM, and mean global signal, the derivatives of these nine regressors, and the squares of these 18 terms. Following confound regression and filtering, the first and last 50 frames of the time series were discarded. Furthermore, a spike regressor was added for each fMRI frame exceeding a motion threshold (0.25-mm root-mean-squared displacement). This confound strategy has been shown to be effective in reducing motion-related artifacts (Parkes et al., 2018). For validation, we also preprocessed the data using aCompCor (Behzadi, Restom, Liau, & Liu, 2007). These data were linearly detrended, band-pass filtered, and trimmed identically to the previous strategy. This confound regression included five high-variance signals estimated from the CSF and white matter each (10 total), as well as six motion estimates, their derivatives, and the squares of these 12 terms. This strategy did not incorporate spike regressors. Following preprocessing and nuisance regression, residual mean BOLD time series at each node was recovered using Connectome Workbench. All data were visually inspected.

### Functional Connectivity and Edge Time Series

Functional connectivity (FC) is generally estimated from fMRI data by computing
the Pearson correlation between the BOLD time series recorded from each node
pair. Hence, each FC estimate represents a linear similarity between the
respective time courses, interpreted as their mutual statistical dependence. It
is, by definition, a nondirected, noncausal metric that does not distinguish
between node pairs that are structurally (anatomically) coupled or uncoupled.
All node pairs maintain nonzero FC, and FC estimates may be negative or
positive. In a system composed of N nodes, the system’s FC matrix has
dimensions [*N* × *N*], due to symmetry
with a total of *K* = (*N*^{2} − *N*)/2 unique entries (all node pairs *i*, *j* with *i* ≠ *j*).

*x*

_{i}is the mean of

*x*(and applied analogously for

*y*) and

*s*

_{x}= $1n\u22121\u2211i=1nxi\u2212X-2$ is the standard deviation of

*x*(and applied analogously for

*y*), and

*n*is the number of observations (for time series data the number of observations is equal to the number of time points). Thus, there are three steps involved in this computation: the conversion of two node time series to z-scores (note that

*z*(

*x*

_{i}) = $xi\u2212X-sx$), forming their product to create a single time series for node pair

*i*,

*j*, and finally forming the mean of this pairwise time series to yield the FC estimate for the node pair (cf. Figure 1). This procedure, when repeated for all pairs of nodes, results in a node-by-node correlation matrix, that is, an estimate of FC. Following an approach developed in prior work (Esfahlani et al., 2020; Faskowitz et al., 2020; Jo et al., 2020), we may omit the final averaging step and retain the time series for each node pair. Since each node pair subtends a unique network edge, we refer to this construct as edge time series. The mean of each edge time series is equal to the corresponding node pair’s FC. Rather than collapsing all time steps into a single scalar FC estimate, omitting the averaging step effectively unwraps FC into a set of edge time series that track co-fluctuation between each node pair. At each point in time, the edge time series report a product of two z-scored BOLD signals; the product is positive if both signals are above or below their respective zero mean, and negative otherwise. The amplitude of their product varies with the joint amplitudes of the two signals. The full set of edge time series comprises a matrix of dimension [

*K*×

*T*], with

*K*equal to the number of unique edges and

*T*equal to the number of time points.

At each moment in time, the amplitude of the co-fluctuations along all edges can
be computed as the root sum square, denoted *RSS*. This metric
takes on high amplitude when edgewise co-fluctuations, on average, are high
(either positively or negatively), and it takes on low amplitudes when
co-fluctuations are low (again, irrespective of their sign). Prior work has
utilized *RSS* to track co-fluctuation amplitudes and stratify or
order time points according to their magnitudes. High-amplitude time points
coincide with intermittent and recurrent patterns of network activity and
connectivity (Esfahlani et al., 2020).

### Bipartitions and Agreement Matrix

Edge time series can be converted to binary form, by applying a threshold at the zero crossings that retains only if co-fluctuations were positive or negative. Positive co-fluctuations occur if and only if two signals both exhibit above mean (positive z-score) amplitudes or if both exhibit below mean (negative z-score) amplitudes. Negative co-fluctuations occur when the two signals deviate in opposite directions. A simple extension of this fact is that, on each time step, positively co-fluctuating node pairs split the network into exactly two communities that are fluctuating negatively with respect to each other. This obligatory two-community split results in a bipartition of the network. The network’s time evolution may be represented as a sequence or set of such bipartitions. Adopting bipartitions largely removes information on signal amplitudes (noting that the bipartition does depend on standardizing the individual node time series) while creating a compact description of the original FC as a set of finely resolved modular partitions, without the need to perform computational community detection. Communities are directly evident from the binary edge time series.

*l*is the number of samples,

*C*

_{k}is the partition of the

*k*-th sample, |

*C*

_{ks}| is the number of nodes in clusters of partition

*C*

_{k}, and

*N*is the number of nodes. In applications to bipartitions from time series,

*l*is equal to the number of time points

*T*. Subtracting the constant null results in agreement matrices that contain negative entries for all node pairs where the observed frequency of co-classification is smaller than that expected under the adopted null model. The subtraction step does not change relative order of frequencies in the agreement matrix and hence has no impact on correlations or similarity metrics computed against FC.

### Bipartition Similarity

*M*and

*M*′ indicate the two partitions,

*m*and

*m*′ indicate modules belonging to the two partitions, and

*P*(

*m*,

*m*′) = $nmm\u2032n$ with

*n*

_{mm′}corresponding to the number of nodes that are members of module

*m*as well as module

*m*′ (Rubinov & Sporns, 2010). In the case of bipartitions, other metrics such as the cosine similarity or the Jaccard distance are also possible. In practice, all these metrics give highly similar results. We adopt the mutual information as the principal metric for assessing similarity between pairs of bipartitions.

### Modularity Maximization

Modularity maximization is a commonly used approach for detecting communities in brain networks (Newman & Girvan, 2004; Sporns & Betzel, 2016) that attempts to partition a network into nonoverlapping communities such that the observed density of connections within subnetworks maximally exceeds what would be expected by chance. The choice of null models should reflect the nature of the network data, which in our case is a correlation matrix (MacMahon & Garlaschelli, 2015). We adopt a constant null (the Potts model; Traag, Van Dooren, & Nesterov, 2011) and retain the full FC matrix, including its negative entries, for the purpose of community detection by applying the Louvain algorithm. Louvain bipartitions are identified by first scanning a wide range of the resolution parameter, selecting upper and lower limits within which a two-community structure appears, followed by a finer sampling of the range to retrieve bipartitions (1,000 samples).

### Gradients

So-called gradients, when computed from FC matrices, represent major connectivity
modes (node vectors) that can be mapped back onto the original node set, for
example the surface of the cerebral cortex. Following a standard workflow (de
Wael et al., 2020), after starting from
an FC matrix, one first derives an affinity matrix that essentially represents a
node-wise distance matrix of size [*N* × *N*]. Here, we compute the affinity matrix, from the full
unthresholded FC, as the cosine similarity for each node pair excluding their
mutual connections. The first principal component of the affinity matrix is
retained for purposes of analysis and comparison. It is identical to the largest
eigenvector of the FC matrix. Other variants for computing the affinity matrix
may include additional steps such as thresholding or alternative distance
transforms. These variants have little to no impact on the relevant findings
reported in this article.

### Individual Differences

*N*×

*N*matrix, from which the upper triangle is converted into a vector of (

*N*

^{2}−

*N*)/2 elements. One such vector is derived for all 95 participants and for each of the four resting-state runs. For each pair of runs (six pairs in total), we compute the Pearson correlation of these vectors across all pairs of participants. The resulting matrix carries similarities between pairs of participants on different runs, with the main diagonal of the matrix representing self-similarity (same participant, different run). The differential identifiability (Amico & Goñi, 2018) is computed as

_{self}corresponding to the mean of the correlations on the main diagonal and $r^$

_{non-self}corresponding to the mean of the off-diagonal correlations. The accuracy is computed as the number of times (out of 95) that the correlation value on the main diagonal (self-similarity) is maximal along the matrix rows and columns. The accuracy estimate is the mean of these two numbers.

Many different criteria can be employed for selecting specific subsets of frames
to compute means of edge time series and agreement matrices, both of size
[*N* × *N*]. Identifiability and
accuracy are computed from such subsets of frames in the same manner as
described above. The selection criterion employed here involves three steps: (a)
computing the mutual information between each frame’s bipartition and all
of the 63 templates in the basis set; (b) normalizing by converting the
resulting 63 *MI* values (per frame) to z-scores; (c) selecting,
for each of the 63 templates, a fraction (here, the top 10%) of frames
with the highest z-scores across a single run. Mean edge time series and
agreement matrices are derived from this subset, representing 10% of the
data. Resulting measures of identifiability and accuracy are compared against a
distribution resulting from several independent realizations (here, 20) of a
null model. In the null model, the timings of the selected frames are shifted by
a random offset (chosen randomly and uniformly between 1 and 1,100). Once the
offset is chosen, frames are rotated in time by the chosen amount. This approach
largely preserves the spacing (autocorrelation) of the frames in the selected
subset.

### Optimization

The purpose of the optimization approach pursued in this study was to discover
sets of bipartitions of *N* nodes, comprising *P* instances, that approximate the observed bipartitions’ (size
[*N* × *T*]) agreement matrix. We
adapted a variant of the Metropolis-Hastings algorithm (Metropolis et al., 1953) to generate samples from the very
large distribution of possible sets of bipartitions of size [*N* × *P*], starting from a random sample and then iteratively
creating variants that are either rejected or accepted, moving into the next
iteration. The rejection or acceptance is governed by simulated annealing
(Kirkpatrick et al., 1983) and by an
objective function *D*, taken here to be a measure of the
distance between the optimized and observed agreement matrix. New variants are
accepted if the objective function improves (lower distance) or if the annealing
criterion *e*^{−ΔD/Temp} > *R*(0, 1) is fulfilled, where *Temp* refers to a simulated “temperature” and *R*(0, 1)
is a random number uniformly drawn from the [0, 1] interval. Essentially, the
annealing criterion allows suboptimal variants to pass, as a function of the
current temperature. The temperature decays exponentially as a function of the
number of iterations, *h*, as *Temp* = *T*_{0}*T*_{exp}^{h}.
Temperature parameters *T*_{0}, *T*_{exp}were selected such that
stable solutions near the global minimum (distance of zero) emerged in
reasonable time. The initial conditions were chosen as sets of completely random
bipartitions, with equal probability (“flipping a coin”) on all
node co-assignments. On each step of the optimization, a single element in a
single bipartition (both chosen at random) was flipped. Note that this
optimization procedure does not implement a true generative process for
bipartitions, as they are not derived from time series data and hence contain no
information on temporal sequences. While more realistic scenarios for
discovering optimally matching sets of bipartitions are conceivable, they were
not pursued in the current study.

Three different formulations of the objective function were tested, all computed
from the agreement matrix’s *K* unique (upper triangle)
elements, with highly reproducible results: (a) the Pearson correlation; (b) the
rank-order correlation (Spearman’s *ρ*); and (c)
the cosine similarity. Optimizations were also carried out by substituting the
observed agreement matrix with the observed FC matrix in the objective function,
with near-identical outcomes.

The number of bipartitions in the optimized set, *P*, is a free
parameter. Different settings of *P*, varying the number of
bipartitions from 1,100 (matching the number of experimental time steps *T*) down to 11, were explored. Smaller values of *P* yield more compact optimized sets, while also resulting
in less accurate matches between the observed and the optimized agreement
matrix.

## ACKNOWLEDGMENTS

Data were provided, in part, by the Human Connectome Project, WU-Minn Consortium (principal investigators: D. Van Essen and K. Ugurbil; 1U54MH091657), funded by the 16 National Institutes of Health (NIH) institutes and centers that support the NIH Blueprint for Neuroscience Research and by the McDonnell Center for Systems Neuroscience at Washington University.

## SUPPORTING INFORMATION

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00182. Matlab code and auxiliary data files demonstrating key concepts introduced in this paper are available at https://www.brainnetworkslab.com/s/bipartitions-code-package.zip (Sporns, Faskowitz, Teixeira, Cutts, & Betzel, 2021a). A movie of edge time series and RSS amplitudes can be downloaded at https://www.brainnetworkslab.com/s/ets_movie_1.mp4 (Sporns, Faskowitz, Teixeira, Cutts, & Betzel, 2021b).

## AUTHOR CONTRIBUTIONS

Olaf Sporns: Conceptualization; Formal analysis; Investigation; Methodology; Software; Visualization; Writing – original draft; Writing – review & editing. Joshua Faskowitz: Conceptualization; Formal analysis; Software; Writing – review & editing. Sofia Teixeira: Conceptualization; Methodology; Writing – review & editing. Sarah A. Cutts: Conceptualization; Formal analysis; Writing – review & editing. Richard Betzel: Conceptualization; Formal analysis; Software; Writing – review & editing.

## FUNDING INFORMATION

Richard Betzel, Indiana University Office of the Vice President for Research, Emerging Area of Research Initiative, Award ID: Learning: Brains, Machines, and Children. Richard Betzel, National Science Foundation (http://dx.doi.org/10.13039/100000001), Award ID: 2023985. Joshua Faskowitz, National Science Foundation (http://dx.doi.org/10.13039/100000001), Award ID: 1342962. Sofia Teixeira, Fundação para a Ciência e a Tecnologia, Award ID: UIDB/50021/2020.

## TECHNICAL TERMS

- Edge time series:
As used here, a time series of co-fluctuations recorded between two nodes (on one edge). In general, any edgewise metric of co-activity, flow, or influence that is measured across time can result in edge time series.

- Bipartition:
A partition of a set of nodes forming a network into exactly two communities. The sizes of the two communities equal the size of the network.

- Co-fluctuation:
As used here, the product of the activation time series of two nodes, converted to standard scores and aligned in time.

- Agreement matrix:
Sometimes also called co-classification or co-assignment matrix. Matrix entries record the frequency, over an ensemble of partitions, with which a node pair is co-assigned to the same network community.

- Template basis set:
Starting from a partition of a network into nonoverlapping communities, a template basis set comprises all possible bipartitions that divide the network into two communities, along the boundaries of the initially chosen partition.

- Identifiability:
A measure that quantifies the relation between similarities in functional connectivity across different imaging sessions for same versus different individual participants.

## REFERENCES

## External Supplements

## Author notes

Handling Editor: Edward Bullmore

DOI:https://doi.org/10.1038/s41598-019-55738-y,PMID:31848397,PMCID:PMC6917755