A mathematical perspective on edge-centric brain functional connectivity — TEI-tagged view

Leonardo NovellipersName, Adeel RazipersName (2022date). “A mathematical perspective on edge-centric brain functional connectivity.” Nature Communicationstitle·j.

Object of study — edge-centric measures of fMRI functional connectivityterm·mth and their relation to the node FCterm·mth

Research question — what information edge-centric measures provide beyond the node functional connectivityterm·mth

Contribution — a mathematical account showing the main edge-centric findings follow from the static node FCterm·mth

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

A mathematical perspective on edge-centric brain functional connectivityterm·mth Abstract Edge time seriesterm·key are increasingly used in brain imaging to study the node functional connectivityterm·mth (nFCterm·mth) dynamics at the finest temporal resolution while avoiding sliding windowsterm·mth. Here, we lay the mathematical foundations for the edge-centric analysis of neuroimaging time series, explaining why a few high-amplitude cofluctuationsterm·dom drive the nFCterm·mth across datasets. Our exposition also constitutes a critique of the existing edge-centric studies, showing that their main findings can be derived from the nFCterm·mth under a static null hypothesis that disregards temporal correlations. Testing the analytic predictions on functional MRI data from the Human Connectome Projectterm·key confirms that the nFCterm·mth can explain most variation in the edge FCterm·mth matrix, the edge communities, the large cofluctuationsterm·dom, and the corresponding spatial patterns. We encourage the use of dynamic measures in future research, which exploit the temporal structure of the edge time seriesterm·key and cannot be replicated by static null modelsterm·key.

Introduction Functional connectivityterm·mth (FCterm·mth) refers to patterns of statistical dependence in brain activity, such as the blood oxygen level-dependent (BOLD) signal measured via functional magnetic resonance imaging (fMRI). Static FCterm·mth is traditionally calculated over the course of an entire scan session and it is an established technique of modern neuroimaging 1 , 2 , with individual differences linked to brain disorders and cognitive states 3 , 4 . On the other hand, time-varying FCterm·mth refers to time-resolved fluctuations in FCterm·mth, typically estimated by fitting dynamic models 5 – 8 or by sliding windowsterm·mth 9 , 10 . Differences in time-varying FCterm·mth at rest are also associated with a wide range of cognitive and behavioural traits, as well as psychiatric and neurological conditions (see ref. 11 for a recent review of this rapidly growing field).

To avoid the temporal blurringterm·dom caused by sliding windowsterm·mth 9 , it is possible to analyse BOLD fluctuations at the resolution of single frames. Coactivation patternsterm·mth (CAPsterm·mth) 12 and point-process analyses 13 are early examples of this approach. Initially using a single seed region, both found that static seed-based FCterm·mth maps can be reliably approximated by averaging only a few high-amplitude frames. Whole-brain extensions of both methods have soon followed 14 , 15 . More recently, edge-centric approaches have generated excitement as they also go beyond single seeds and additionally decompose the entire FCterm·mth matrix into its frame-wise contributions by omitting traditional time averaging over the length of the experiment 16 , 17 . Such temporal unwrapping of the FCterm·mth results in a large number of edge time seriesterm·key, each capturing the moment-to-moment cofluctuationsterm·key of a pair of brain regions. Going one step further, one can measure the similarity between every pair of edges to obtain a large matrix, referred to as edge FCterm·mth (eFCterm·mth), as opposed to the traditional node-centric FCterm·mth (nFCterm·mth). The eFCterm·mth is shown to be replicable, stable within individuals across multiple scanning sessions and reliable across datasets 17 . Furthermore, clustering the eFCterm·mth yields overlapping brain communities that could better suit the study of aspects of cognition and behaviour that transcend traditional disjoint brain parcellations. While CAPsterm·mth focus on patterns of BOLD activity, edge-centric approaches focus on patterns of cofluctuationsterm·key between all regional pairs. However, the analogous finding is that the static nFCterm·mth can be faithfully approximated by averaging a few instantaneous connectivity patterns. These are characterised by simultaneous large cofluctuationsterm·dom across all node pairs, measured as the root-sum-of-squaresterm·mth (RSSterm·mth) over all the edge time seriesterm·key values corresponding to the same frame. Such brief, intermittent, and high-amplitude cofluctuationsterm·dom drive the nFCterm·mth and the network structure over these time points contributes disproportionately to the overall modularity of the functional brain network 16 .

Although the edge-centric decomposition of nFCterm·mth into its frame-wise constituents is mathematically exact, a comprehensive treatment of the statistical properties of the edge-centric measures is lacking and there is a consensus on the need for appropriate null modelsterm·key 17 , 18 . A rigorous mathematical study is especially important since several widely-acknowledged publications have warned about the dangers of extracting structure from noise when studying static or time-varying FCterm·mth, often using minimal null modelsterm·key to reproduce existing results 19 – 24 . The warnings concerning sampling variability are particularly relevant to edge-centric methods as they represent an extreme case of a single-frame sliding-window approach. High-amplitude cofluctuationsterm·dom can be observed in temporally-uncorrelated synthetic time series such that accounting only for static spatial correlations is sufficient to replicate key empirical findings 16 (Fig. S4) . This observation has been interpreted as further evidence that large cofluctuationsterm·dom are not fMRI artefacts. However, it arguably raises an equally pressing conceptual concern: what information do current edge-centric measures provide beyond the nFCterm·mth, if any?

Here, we tackle this question mathematically and present a theoretical explanation for the widespread occurrence of large cofluctuationsterm·dom across datasets and why a few large events drive the nFCterm·mth. This explanation rests on fundamental properties of subexponential distributions 25 . Further mathematical derivations clarify how the nFCterm·mth eigenvaluesterm·key shape the RSSterm·mth distribution and how the leading nFCterm·mth eigenvectors underpin the spatial correlation patterns expressed during high-amplitude events. The influence of functional modules on the eigenvalueterm·key distribution could explain why these events disappear when the modular structure is disrupted, as recently reported in ref. 26 . Finally, we analytically show that the eFCterm·mth matrix, the edge communities, the large cofluctuationsterm·dom, and the corresponding brain activity modes can all be predicted from the nFCterm·mth without recourse to the edge-centric formulation. Many of these derivations are based on the null hypothesis of i.i.d. Gaussian variables that only takes into account the observed (static) spatial correlations and ignores temporal features. Under this assumption, and invoking results from random matrix theory 27 , the edge time seriesterm·key variability is described by the sampling distribution of the nFCterm·mth, known as the Wishart distributionterm·key 28 . Testing the analytic predictions using fMRI data from the Human Connectome Projectterm·key (HCP) 29 shows that the null modelterm·key is sufficient to replicate the vast majority of existing edge-centric features both qualitatively and quantitatively, as well as foundational properties of CAPsterm·mth.

Results We present six main results showing that the existing findings based on edge time seriesterm·key 16 , 17 can be derived from the static nFCterm·mth under the null hypothesis of i.i.d. multivariate Gaussian variables that preserve the observed static spatial correlations but not the temporal ones. The theoretical predictions were empirically tested using a HCP dataset comprising 100 subjects, preprocessed with the current standard HCP pipeline, both with and without global signal regressionterm·key (GSR). The detailed derivations of the presented equations are in the ‘Methods’ section, reserving the current section for a concise account of the key results.

The edge FCterm·mth matrix can be derived analytically from the node FCterm·mth The eFCterm·mth was introduced in ref. 17 to quantify linear interactions between edges. For each pair of brain regions, an edge time seriesterm·key is computed as the element-wise product of the two regional z-scored time series. Thus, the values of each edge time seriesterm·key represent the instantaneous cofluctuationterm·key magnitudes between the corresponding pair of brain regions. The eFCterm·mth is the edge-by-edge matrix obtained by computing the inner products between all pairs of edge time seriesterm·key, normalised to the [−1, 1] interval.

Our first result is that the eFCterm·mth can be analytically derived from the nFCterm·mth under the static Gaussian null hypothesis (Fig. 1 a, b). As shown in Eq. ( 11 ) of the ‘Methods’ section, the ( j k , l m ) entry of the eFCterm·mth matrix is obtained as a sum of pairwise products between the ( j , k , l , m ) entries of the nFCterm·mth matrix, divided by a normalisation factor: 1 eFCterm·mth j k , l m = r j k r l m + r j l r k m + r j m r k l 1 + 2 r j k 2 1 + 2 r l m 2 . Using the HCP dataset, the predicted eFCterm·mth achieves an average Pearson correlation of r = 0.93 with the empirical eFCterm·mth ( r = 0.88 if GSR is applied). The distributions across 100 unrelated HCP subjects are shown in Fig. 1 c. This is a significant improvement on the linear regression approach adopted in ref. 17 , which achieved an average Pearson correlation of r = 0.72 on pairs of edges not sharing any nodes, but performed poorly otherwise ( r = 0.06). Moreover, by revealing the mathematical relationship between eFCterm·mth and nFCterm·mth, Eq. ( 1 ) explains why the eFCterm·mth is highly replicable, stable within individuals across multiple scan sessions and consistent across datasets—it will be as long as the nFCterm·mth is. Fig. 1 The edge functional connectivityterm·key (eFCterm·mth) can be derived analytically from the node functional connectivityterm·mth (nFCterm·mth) under the static Gaussian null modelterm·key. a The eFCterm·mth of the first subject from the Human Connectome Projectterm·key (HCP) fMRI dataset is plotted as an example. This matrix scales with the fourth power of the parcellation size: in the current case of 200 parcels, it has ~400 million elements. The ( j k , l m ) entry of the eFCterm·mth matrix (black square at the intersection of the black lines) is computed as the inner product of two edge time seriesterm·key j k and l m , normalised to the [−1, 1] interval. b Under the static null modelterm·key, the eFCterm·mth of each subject can be derived analytically as a sum of pairwise products between the ( j , k , l , m ) entries of the corresponding nFCterm·mth matrix as per Eq. ( 1 ). Since the nFCterm·mth is symmetric, only six entries are necessary to represent all pairwise combinations (marked by each of the black squares at the intersections of the black lines). c Distribution of Pearson correlation coefficients between the empirical and null eFCterm·mth computed over 100 HCP subjects. Results are reported with and without global signal regressionterm·key (GSR) and the mean correlation values are indicated by the dashed lines.

The edge communities can be predicted from the nFCterm·mth Clustering the eFCterm·mth via the k-means algorithm was used to identify communities of co-fluctuating edges. These were then mapped back to individual nodes to obtain overlapping regional communities, offering a new way to study aspects of cognition and behaviour that transcend traditional disjoint brain parcellations 17 . Given the high similarity between the empirical and null eFCterm·mth demonstrated in the previous section, it would be reasonable to expect that the edge communities be well recovered by the null modelterm·key by applying the same community detection algorithm to the predicted eFCterm·mth. Indeed, the agreement between the empirical and predicted community assignments is exact on 84% of the 19,900 edges, and 74% if GSR is applied (Fig. 2 a, b). At the node level, the similarity is even stronger. In ref. 17 , the similarity between two nodes, referred to as “edge cluster similarity”, is measured as the fraction of all edge pairs starting from those two nodes (and reaching the same target) that are clustered together. Using the same procedure, the predicted edge cluster similarity achieves a Pearson correlation coefficient of r = 0.96 with the empirical one, and r = 0.95 if GSR is applied (Fig. 2 c, d). Note that the rows and columns of the matrices in Fig. 2 have been rearranged to match the ordering of the 16 networks used in ref. 17 (Fig. 6 ) to facilitate a visual comparison. Given that the analysis is performed on the same HCP dataset, the small differences in the empirical results are most likely due to the different preprocessing pipelines (here we use the preprocessed data made available by the HCP; see the ‘Methods’ section for details). This adds to the evidence that edge-centric measures are not strongly dependent on the preprocessing approach, and similarly speaks to the robustness of the null modelterm·key predictions. Fig. 2 The edge communities can be predicted from the nFCterm·mth under the static Gaussian null modelterm·key. a Empirical edge communities obtained by k-means clustering of the eFCterm·mth computed on the HCP dataset. The edge labels are reshaped into an N x N matrix where each ( i , j ) entry represents the community label of the edge linking i to j . For each edge, the representative community assignment shown here is the statistical mode (i.e., the most common label) across 100 HCP subjects. b Predicted edge communities obtained by k-means clustering of the eFCterm·mth predicted by the static null modelterm·key. The agreement with the empirical labels in panel a is exact on 84% of all 19900 edges (74% if GSR is applied). c Empirical edge cluster similarity on the HCP dataset. The similarity of edge communities involving nodes i and j is computed as the fraction of matching elements between the corresponding two rows of the edge community matrix in panel a . d Predicted edge cluster similarity based on the static null modelterm·key. The Pearson correlation with the empirical edge cluster similarity in panel c is r = 0.96 ( r = 0.95 if GSR is applied).

So far, the results in this subsection have been based on simulations. Let us now see how a mathematical derivation can provide further insight into these edge and node similarities. First, recall that the edge communities are obtained by clustering the eFCterm·mth (e.g. via k-means). The main obstacle to a full analytic approach is that the outcome of stochastic clustering algorithms cannot be entirely predicted from their input; however, since they are usually based on a distance metric, it is reasonable to expect that the smaller the distance between two rows of the eFCterm·mth, the higher the probability that the corresponding edges would be clustered together. Accordingly, we define the distance between two edges ( j k ) and ( j ′ k ′ ) as the ℓ 1 norm of the difference between the corresponding rows of the eFCterm·mth. As shown in Eq. ( 12 ) of the ‘Methods’ section, this edge distance simplifies to 2 d j k , j ′ k ′ = ∑ l , m = 1 N 3 z j ( z l ⊤ z m ) z k ⊤ − z j ′ ( z l ⊤ z m ) z k ′ ⊤ , where, for a brain region i , the row vector z i is its z-scored BOLD signal and z i ⊤ denotes its transpose. What does this imply for the node communities? The similarity between two nodes i and j was measured in ref. 17 as the fraction of all edge pairs starting from i and j (and reaching the same target) that are assigned to the same cluster. Here, the analogous analytic step is to compute the distance between nodes i and j as the sum of the distances between the edges starting from them. Crucially, the resulting node distance can be expressed in terms of BOLD signal correlations (see Eq. ( 13 ) in the ‘Methods’ section): 3 d i , j ≤ c 1 − r i j 1 2 , where c denotes a constant term, independent of i and j . This implies that the edge-cluster similarity 17 between nodes i and j can be predicted from the corresponding r i j entry of the nFCterm·mth, avoiding the memory-intensive computation of the eFCterm·mth and computationally-intensive clustering algorithms (the space complexity of the eFCterm·mth is O ( N 4 ) , i.e., it scales with the fourth power of the number of regions and requires over a terabyte of memory for fine brain parcellations—for each subject). Using the HCP dataset to test this prediction shows that the nFCterm·mth alone achieves an average Pearson correlation of r = 0.76 with the empirical edge cluster similarity matrix shown in Fig. 2 c. Once again, the static, node-centric, second-order features of the BOLD signal are sufficient to replicate key findings that appear at first to rely on the specific temporal sequence of BOLD cofluctuationsterm·key at the single-frame resolution.

The null modelterm·key reproduces the high similarity of the top RSSterm·mth frames to the nFCterm·mth Let us now consider the root-sum-of-squaresterm·mth (RSSterm·mth) of the edge time seriesterm·key introduced in ref. 16 . The RSSterm·mth is a univariate time series defined as the Euclidean norm of the edge time seriesterm·key vector at each frame. In other words, the RSSterm·mth peaks when all the cofluctuationsterm·key (i.e., the edge time seriesterm·key) are simultaneously high in absolute value, either positive or negative. The key finding in ref. 16 is that only a small fraction of frames exhibiting large RSSterm·mth values are required to explain a significant fraction of variance in the nFCterm·mth, as well as the network’s modular structure 30 . Both results are perfectly reproduced by the static null modelterm·key, as shown in Fig. 3 a, b and Supplementary Fig. 1 . What is particularly remarkable is that the timing of the high-amplitude RSSterm·mth events produced by the null modelterm·key are arbitrary, and yet a small fraction of frames corresponding to these large cofluctuationsterm·dom is still sufficient to explain the observed nFCterm·mth. Furthermore, Fig. 3 c shows that the null modelterm·key frames with the largest RSSterm·mth also exhibit high similarity to the empirical frames with the largest RSSterm·mth from the HCP dataset—occurring at entirely different times. Note that, unlike temporal dependencies, spatial correlations are necessary to reproduce the results: if the null modelterm·key is chosen to be both temporally and spatially uncorrelated, high-RSS frames are no more similar to the empirical nFCterm·mth than low-RSS frames (see Supplementary Fig. 2) . Another interesting note is that high-RSS frames exhibit the strongest average correlation with all other BOLD frames and this average similarity decreases with the RSSterm·mth magnitude, both in empirical and null cases (see Supplementary Fig. 3) . Fig. 3 The static Gaussian null modelterm·key can reproduce the strong correlation between high-amplitude cofluctuationterm·dom patterns and the nFCterm·mth. a Only a small fraction of frames exhibiting the largest cofluctuationsterm·key root-sum-of-squaresterm·mth (RSSterm·mth) are required to explain most of the nFCterm·mth variance. The similarity is computed as the Pearson correlation coefficient between the nFCterm·mth and the average FCterm·mth estimated from the top and bottom 5% of the total frames. Each point corresponds to one of 100 unrelated subjects from the Human Connectome Projectterm·key (HCP) dataset, with boxes indicating the quartiles and whiskers length specified as 1.5 times the interquartile range. b The same results hold more generally when the frames are ordered according to the corresponding RSSterm·mth amplitude, either in descending or ascending order. Here, the curves represent the average similarity over 100 subjects. c The findings do not depend on the timing of the high-amplitude RSSterm·mth events: the top frames generated by the null modelterm·key exhibit high similarity to the top frames of the real HCP data, which occur at different times.

A theoretical explanation for these findings will be provided in the ‘nFCterm·mth eigenvectors underpin spatial patterns of high BOLD activity’ section and supported by detailed derivations in the ‘Methods’ section. Interestingly, most of these points were also reported in ref. 16 (Fig. S4) , where they were taken as evidence that large RSSterm·mth events are not fMRI artefacts. While settling that methodological issue, these observations raise a conceptual concern: if matching the timing of the RSSterm·mth events is not essential and the results can be replicated by the static null modelterm·key, does the edge-centric approach provide any information about the time-varying connectivity that cannot be explained by the static nFCterm·mth? We will address this question in the next section, by examining the statistical properties of the RSSterm·mth.

The RSSterm·mth distribution is determined by the nFCterm·mth eigenvaluesterm·key Having established that the large RSSterm·mth events are not an exclusive feature of neural signals, let us investigate how their ubiquitous appearance across datasets can be analytically explained and why the corresponding frames account for the largest fraction of variance in the nFCterm·mth. As a first step, the RSSterm·mth can be computed as the squared Euclidean norm of the z-scored BOLD signal z (full derivations provided in the ‘Methods’ section, see Eqs. ( 14 ) to ( 16 ): 4 RSSterm·mth all ( t ) = z ( t ) 2 . The intuition behind this equivalence is that summing over the pairwise products of the elements of a vector is the same operation that is performed when squaring a polynomial; in this case, the vector is the squared BOLD signal, its pairwise products are the squared edge time seriesterm·key, and the polynomial is the squared Euclidean norm. The key message is that, although it was introduced in ref. 16 as the Euclidean norm of the edge time seriesterm·key, the RSSterm·mth is mathematically equivalent to the squared Euclidean norm of the BOLD signal, i.e., a measure of the overall BOLD signal amplitude at each time step, akin to the variance. We can then proceed without resorting to the edge time seriesterm·key, which is not only convenient in practice but also shifts the conceptual focus back to the BOLD time series—which are more readily interpretable. For a large family of common (sub-Gaussian) distributions, the squared Euclidean norm of a random variable (RV) is heavy-tailed (more specifically, it is subexponential 31 ). The RSSterm·mth being a squared norm, large cofluctuationsterm·dom are then to be expected due to its heavy-tailed distribution (see Eq. ( 28 )), offering an explanation for the large RSSterm·mth peaks observed in the BOLD time series.

In the specific case of Gaussian variables (i.e., the null hypothesis), the RSSterm·mth can be expressed as a sum of N independent Gamma ( k = 1 / 2 , θ = 2 λ i ) variables, each related to an eigenvalueterm·key of the nFCterm·mth matrix ( λ i , i = 1…, N ). This is summarised by the moment-generating function in Eq. ( 27 ). The largest eigenvaluesterm·key capture the distribution tail and including smaller eigenvaluesterm·key provides an increasingly complete characterisation of the empirical RSSterm·mth distribution (Fig. 4 b). This distribution can be used for testing the statistical significance of the empirical RSSterm·mth observed in the HCP dataset against the static null hypothesis of spatially correlated noise. Figure 4 a illustrates the convergence of the empirical RSSterm·mth distribution to the null distribution as more and more time frames are observed (i.e., over longer fMRI sessions). When all the 1200 time frames available in the HCP data are utilised, the null hypothesis cannot be rejected for 58% of the participants at a 5% significance level and on 90% of the participants after Bonferroni correction for multiple comparisons ( p -values given by the two-sided Kolmogorov–Smirnov test). Fig. 4 Properties and convergence of the RSSterm·mth distribution derived analytically under the static null hypothesis. a The empirical RSSterm·mth distribution increasingly approximates the null RSSterm·mth distribution as more realisations (frames) are observed. On the horizontal axis, the cumulative distribution function (CDF) of the null RSSterm·mth is evaluated for each empirical RSSterm·mth observation; the CDF of the resulting distribution is plotted on the vertical axis. Identical distributions would produce points on the diagonal (dotted line). b The nFCterm·mth eigenvaluesterm·key shape the null RSSterm·mth distribution. The null RSSterm·mth distribution increasingly approximates the empirical RSSterm·mth distribution as more eigenvaluesterm·key are considered, with the largest eigenvaluesterm·key explaining most of the variance. On the horizontal axis, the cumulative distribution function (CDF) of the empirical RSSterm·mth is evaluated for each null RSSterm·mth observation; the CDF of the resulting distribution is then plotted on the vertical axis. Identical distributions would produce points on the diagonal (dotted line).

nFCterm·mth eigenvectors underpin spatial patterns of high BOLD activity High-amplitude BOLD cofluctuationsterm·key were observed to be underpinned by a particular spatial mode of brain activity in which default mode and control networks are anticorrelated with sensorimotor and attentional systems 16 . This particular spatial mode was defined as the first principal component of the BOLD activity and, as such, it can be obtained as the largest eigenvector of the static nFCterm·mth matrix by mathematical equivalence and without recourse to null modelsterm·key (Fig. 5 a). Fig. 5 The static Gaussian null modelterm·key can explain the spatial patterns of BOLD activity and the cofluctuationterm·key patterns that characterise high-amplitude frames. a The spatial mode underpinning high-amplitude cofluctuationsterm·dom is captured by the leading eigenvector of the nFCterm·mth matrix. The 200 brain regions were partitioned into the same 16 networks used in ref. 16 (Fig. 2 e) to facilitate a visual comparison. In addition to a scatter plot of the data, the violin plots show the probability densities and the box plots indicate the quartiles, with the maximum whisker length specified as 1.5 times the interquartile range. The sample sizes of the box plots are the sizes of the 16 networks, that is, n = [16, 15, 6, 14, 17, 6, 12, 10, 14, 16, 10, 19, 15, 6, 12, 12]. b Frames corresponding to large RSSterm·mth values exhibit high similarity to the leading nFCterm·mth eigenvectors. The similarity is measured as the Pearson correlation between the instantaneous FCterm·mth estimate from a single frame and each of the estimates from the four leading nFCterm·mth eigenvectors. c Higher principal component (PC1) coefficients are associated with large RSSterm·mth events. In addition to a scatter plot of the data, the box plots indicate the quartiles and the whisker length is specified as 1.5 times the interquartile range (sample size n = 6000). Compare this figure with the results published in ref. 16 (Fig. 2 c).

The only question left to answer is whether the RSSterm·mth is predicted to peak when the BOLD activity aligns with the largest eigenvector. This can be proven to be true even beyond the static null hypothesis (see Eqs. ( 23 ) to ( 24 ) in the ‘Methods’ section). An intuitive understanding can be gained once the RSSterm·mth is seen as the fluctuation of the BOLD signal amplitude over time: high-amplitude frames have a larger variance, which is captured by a larger coefficient of the first principal component since the latter is the vector that aligns with the direction of maximum variance (Fig. 5 c). We can then refine our theoretical understanding of the RSSterm·mth peaks: not only do they occur when the Euclidean norm of the BOLD signal vector is large but, most likely, when it is well aligned with the leading eigenvector of the nFCterm·mth. However, the alignment with the leading eigenvector need not be perfect: in general, large RSSterm·mth values can be expected whenever the BOLD vector is a mixture of the top eigenvectors (Fig. 5 b and Supplementary Fig. 4) . Additional principal components are included as more large-RSS frames are averaged, suggesting why the top 5% frames alone are sufficient for an almost perfect reconstruction of the nFCterm·mth (Fig. 3) . We have thus explained why the nFCterm·mth estimates corresponding to frames with the largest RSSterm·mth exhibit the highest similarity with the nFCterm·mth. Moreover, since the nFCterm·mth features multiple communities, high-RSS frames naturally reflect this property by exhibiting high modularity, as well as higher values than low-RSS frames, which are less similar to the nFCterm·mth (see Supplementary Fig. 1) .

One could speculate that the cofluctuationterm·key patterns corresponding to large nFCterm·mth eigenvectors would be closely related to empirical cofluctuationterm·key patterns obtained by clustering high-RSS frames as in ref. 18 . For example, the leading nFCterm·mth eigenvector has the highest similarity to all BOLD frames and the corresponding cofluctuationterm·key pattern would resemble the most frequently occurring cluster.

A null modelterm·key for binary edge time seriesterm·key It was recently observed that thresholding the edge time seriesterm·key retains most of the information about the nFCterm·mth since averaging the binarised edge time seriesterm·key over time still yields a very accurate approximation of the nFCterm·mth both at the voxel 15 and parcel level 32 . As a premise, the latter finding is empirically replicated here using the HCP dataset, with a resulting Pearson correlation r = 0.98 between the average binarised edge time seriesterm·key and the nFCterm·mth (average correlation over 100 unrelated subjects). How to explain this almost perfect correlation? Furthermore, it was noted that the binary edge time seriesterm·key are highly constrained by the nFCterm·mth 32 . What is the nature of such constraints? Both questions can be answered mathematically under the static Gaussian null hypothesis. Consider two brain regions or parcels j and k with z-scored BOLD activity Z j ( t ) and Z k ( t ) and their corresponding edge time seriesterm·key C j k ( t ) = Z j ( t ) Z k ( t ). The thresholded edge time seriesterm·key C j k ( t ) ¯ is equal to 1 when the cofluctuationsterm·key between j and k are positive, i.e., when Z j ( t ) and Z k ( t ) have the same sign. Under the null hypothesis, Z j ( t ) and Z k ( t ) are normal Gaussian RVs with correlation coefficient denoted as r j k . If the two parcels are uncorrelated ( r j k = 0), their sign is positive or negative with equal probability and can be described mathematically as a Bernoulli(1/2) RV, which is a formal way to model a fair coin draw. The corresponding binary edge time seriesterm·key C j k ( t ) ¯ can also be modelled as a fair coin draw since the two signs are expected to agree half of the times. The situation clearly gets more complicated in the case of 200 correlated parcels as in the HCP dataset considered here. Intuitively, we would expect C j k ( t ) ¯ = 1 if the two parcels were perfectly correlated and C j k ( t ) ¯ = 0 if they were perfectly anticorrelated, i.e., in consistent disagreement. This intuition can be formalised mathematically and extended to all intermediate cases as shown in Eq. ( 30 ) of the ‘Methods’ section. A trigonometric argument proves that C j k ( t ) ¯ can be modelled as a biased coin, i.e., a Bernoulli( p j k ) RV with success probability 5 p j k = 1 2 + arcsin ( r j k ) π . This result reveals the exact analytic relationship between the binary edge time seriesterm·key and the nFCterm·mth, answering the second question. However, it is straightforward to show that Eq. ( 5 ) also predicts the approximated nFCterm·mth obtained by averaging the binary edge time seriesterm·key. This follows from the fact that the time average of ergodic processes (including the null modelterm·key) converges to their expectation, and that the expectation of a Bernoulli RV is equal to its success probability, p j k . Testing this analytic prediction on the HCP dataset confirms its accuracy (see Supplementary Fig. 5 b). Finally, we can explain the strong correlation between the original and approximated nFCterm·mth by noting that p j k in Eq. ( 5 ) is very close to r j k for small values of the correlation (see Supplementary Fig. 5 a), which are the most frequent ones.

Relationship with coactivation patternsterm·mth (CAPsterm·mth) CAPsterm·mth 12 , 33 preceded edge-centric approaches in the study of spatial patterns of BOLD activity at the single-frame resolution. These patterns differ from established resting-state functional networks and it was speculated that they may originate from a neuronal avalanching phenomenon. Both CAPsterm·mth and edge-centric studies report a strong similarity between high-amplitude frames and the nFCterm·mth, and both employ k-means to assign these frames to different clusters 12 , 18 . However, there are important differences between the two methods. CAPsterm·mth are patterns of node coactivations (i.e., they are defined in node space), while the frame-wise FCterm·mth patterns studied in edge-centric analyses are defined in edge space. Moreover, CAPsterm·mth are seed-based, i.e., frames to be clustered are selected based on high activity levels of a single node rather than simultaneous high activity across all nodes. For example, choosing the posterior cingulate cortex (PCC) as the seed, Liu and Duyn 12 observed a strong correlation between the BOLD activity at times where PCC is highly active and the corresponding seed-based FCterm·mth map.

The static Gaussian null modelterm·key employed in this study can predict this strong correlation, which is the fundamental conceptual and practical property guiding the selection of frames that are subsequently clustered into CAPsterm·mth. As shown in Eq. ( 32 ) of the ‘Methods’ section, frames selected based on the high activity of a seed node are expected to exhibit BOLD activity patterns that reflect the corresponding column of the nFCterm·mth. Furthermore, the correlation is directly proportional to the activity level and thus peaks at frames with the highest activity of the seed. This relationship is illustrated in (Fig. 6 a), where frames are sorted in descending order based on the seed node activity. It is important to note that parcels are used instead of voxels, in alignment with the rest of the simulations herein; however, the mathematical results hold true in both cases. What happens when more frames are averaged? The answer is shown in Fig. 6 b and follows intuitively from Fig. 6 a. Let us again consider the frames sorted in descending order; as the activity threshold is lowered (i.e., set to a lower percentile), more frames are averaged. The first frames are best aligned with the seed FCterm·mth vector (i.e., the nFCterm·mth column corresponding to the seed node) and their average quickly converges to the same direction. The alignment (or similarity) reaches a plateau as the middle frames are added to the average since they have low or zero correlation with the seed FCterm·mth vector. Finally, as the last and most negatively-correlated frames are added to the average, the similarity drops sharply. This explains and replicates the first findings of the seminal CAPsterm·mth paper by Liu and Duyn 12 . Fig. 6 The static Gaussian null modelterm·key can explain the spatial similarity between high-amplitude frames based on a seed node and the corresponding correlation map—a core feature of coactivation patternsterm·mth (CAPsterm·mth). a For each seed node (i.e. parcel), the frames are sorted in descending order based on the BOLD activity of the seed (horizontal axis). The Pearson correlation is computed between each frame and the seed-based correlation map (i.e., the column of the nFCterm·mth matrix corresponding to the seed). The curves represent the average similarity over all 200 seed nodes and over 100 HCP subjects. As predicted analytically, the similarity is proportional to the BOLD activity and it decreases as lower-amplitude frames are considered (lower percentiles). The null modelterm·key is in excellent agreement with the empirical results. b Same as in panel ( a ), but the BOLD activity of all the frames above a given threshold (percentile) are averaged before computing the correlation with the seed-based FCterm·mth map. Only a small fraction of high-amplitude frames is required to explain most of the nFCterm·mth variance, reproducing a well-known result in the CAPsterm·mth literature 12 .

To summarise, analytic predictions in the case of CAPsterm·mth differ from the edge-centric case in that the high-amplitude frames are best explained by individual nFCterm·mth columns rather than its eigenvectors. However, there is an intuitive link between the two: if we select frames where a specific node is highly active (CAPsterm·mth approach), the BOLD signal will most likely resemble the corresponding nFCterm·mth column; if we select frames where all nodes are highly active (edge-centric approach), the BOLD signal will most likely resemble all the nFCterm·mth columns, i.e., it will align with the leading nFCterm·mth eigenvector.

Discussion We have presented mathematical proofs and numerical analyses of real HCP data supporting our claim that the static nFCterm·mth is sufficient to replicate the main resting-state edge-centric findings in refs. 16 , 17 both qualitatively and quantitatively, without relying on the edge time seriesterm·key nor any temporal correlations. Specifically, the eFCterm·mth, the edge communities, the edge time seriesterm·key norm (RSSterm·mth) distribution, and the spatial BOLD patterns underpinning large cofluctuationsterm·dom can all be predicted from the nFCterm·mth under the null hypothesis of i.i.d. multivariate Gaussian variables. As further shown, key properties of binary edge time seriesterm·key 32 and CAPsterm·mth 12 can similarly be predicted. The inability to reject the null hypothesis on most of the HCP 100 unrelated subjects does not support the conclusion that these edge-centric metrics provide additional information beyond the nFCterm·mth. These results are not an attempt to disprove the existence of finely timed neural events—they just warn that the evidence provided by fMRI data may not be sufficient to reject simpler explanations, based on the established nFCterm·mth literature, for the edge-centric features studied in refs. 16 , 17 . In fact, previous influential studies have raised similar warnings in the context of sliding-window approaches to time-varying FCterm·mth 19 , 21 , 22 .

However, it would be premature to conclude that the nascent edge-centric approach has no merit, and we acknowledge the fast progress in its development and applications at the time of writing 18 , 26 , 32 . Particularly interesting is the influence of structural modules on the edge cofluctuationsterm·key 26 , which we briefly address in the ‘Methods’ section. The size of the functional modules shapes the spectrum of the nFCterm·mth: larger functional modules allow for larger eigenvaluesterm·key, which underpin the high-amplitude cofluctuationsterm·dom. This offers a mathematical insight into the relationship between modular structure and large cofluctuationsterm·dom, and why the latter disappear if the modular structure is disrupted. While fully addressing the latest (and increasingly large number of) preprints is beyond the scope of this work, it is possible that model-based approaches will reveal the role of edge-centric properties in bridging brain structure and function. Indeed, temporally-unfolded (or point-wise) dependence measures have been instrumental in studying the structure-function relationship in canonical complex systems 34 , 35 ; seeing the edge time seriesterm·key as point-wise mutual information under the Gaussian assumption could create new links to the existing literature.

It would also be unreasonable to assume that the null hypothesis of i.i.d. variables is a good description of the BOLD signal, which is slowly-varying and highly autocorrelated. Thus, the fact that such null modelterm·key is able to replicate the edge-centric features in refs. 16 , 17 is an indication that the temporal structure of the edge time seriesterm·key has not been fully exploited. Indeed, besides the notable cases of the synchronisation of the cofluctuationsterm·key across subjects watching the same movies and the peak-to-peak interval distribution (see Fig. 6) , most of the proposed edge-centric metrics are unaffected by time shuffling (i.e., they are static measures) and can thus be replicated by the i.i.d. null modelterm·key. Our mathematical derivations do not directly apply to sliding-window approaches since windowed correlations are changed by time shuffling, but it is important to remember that many common time-varying FCterm·mth analysis pipelines have intermediate steps that alternately leverage and neglect temporal ordering 11 . For example, one might estimate sliding-window correlations (dynamic stage), apply k-means clustering to the resulting time-resolved FCterm·mth matrices (static stage, since k-means ignores the temporal ordering of the windows), and then evaluate state properties such as dwell times and transition probabilities (dynamic stage). While dynamic measures cannot be predicted from the static nFCterm·mth, it is not unlikely that static null modelsterm·key could reproduce some of the static measures involved, e.g. the state patterns found via k-means, but not their transition probabilities.

Let us note, however, that the role of null modelsterm·key in time-varying FCterm·mth is a matter of current debate 19 , 36 , and not all features that can be explained by null modelsterm·key are clinically irrelevant or to be dismissed. For example, using a small fraction of high-amplitude frames to approximate the nFCterm·mth has been suggested as a way to compress the BOLD signal and alleviate the computational burden of analysing large fMRI datasets without compromising the prediction accuracy 15 . Future studies may find other useful criteria for filtering the frames using the edge time seriesterm·key. As a final contribution, we have analytically shown that the predominance of a few frames in shaping the nFCterm·mth is to be expected: the nFCterm·mth captures the BOLD signal variance, which is a second-order statistic and heavy-tailed (even if the BOLD signal were Gaussian). Therefore, the nFCterm·mth is necessarily shaped by a few tail events corresponding to large amplitude frames. These frames determine the direction of maximum variance and thus the first principal components of the BOLD signal, forming the leading eigenvectors of the nFCterm·mth.

In conclusion, we have laid out the mathematical foundations for the edge-centric FCterm·mth analysis with the goal of informing future studies, in an interplay with empirical observations and simulations. Future work could leverage this mathematical framework and focus on dynamic measures that cannot be easily explained by minimal null modelsterm·key like the one presented herein.

Methods Definition of edge-centric FCterm·mth Functional connectivityterm·mth is defined as the magnitude of the statistical dependence between pairs of brain parcels 37 . This dependence is typically estimated from their time series (here, the BOLD signal) using the Pearson correlation coefficient. Let N be the number of parcels, T be the number of recorded frames, and x i = [ x i (1), …, x i ( T )] be the time series recorded from parcel i , with 1 ≤ i ≤ N . The correlation between two parcels i and j can be computed as r i j = 1 T − 1 ∑ t z i ( t ) z j ( t ) , where z i and z j are their z-scored time series row vectors, i.e., z i = x i − μ i σ i (with μ i and σ i indicating the time-averaged mean and standard deviation). Repeating this procedure for all pairs of parcels results in a node-by-node ( N × N ) correlation matrix R = [ r i j ], which is an estimate of the (node-centric) functional connectivityterm·mth.

The edge time seriesterm·key between two parcels i and j is the vector resulting from the element-wise product of z i and z j , which encodes the magnitude of their cofluctuationsterm·key over time: 6 c i j ( t ) : = z i ( t ) z j ( t ) . We will denote the random variable (RV) associated with the edge time seriesterm·key c i j ( t ) with a capital letter, i.e., C i j ( t ). The column vector of all the N 2 edge time seriesterm·key values at a given time t can be reshaped into a N × N matrix that is an instantaneous estimate of the dynamic functional connectivityterm·mth based on a single frame. This matrix is symmetric since each edge time seriesterm·key is independent of the order of the node pair it involves, i.e., c i j ( t ) = c j i ( t ), ∀ i , j . For this reason, only the upper-triangular portion is typically computed, for a total of N ( N − 1)/2 edge time seriesterm·key instead of N 2 .

It is also possible to go one step further and estimate the statistical dependence between each pair of edge time seriesterm·key, where each edge corresponds to a pair of parcels. This fourth-order statistic results in a large N ( N − 1 ) 2 × N ( N − 1 ) 2 matrix named edge functional connectivityterm·key (eFCterm·mth) 17 , with normalised entries defined as 7 eFCterm·mth j k , l m : = ∑ t c j k ( t ) c l m ( t ) ∑ t c j k ( t ) 2 ∑ t c l m ( t ) 2 .

Null hypothesis Let z be the N × T (parcels × frames) matrix of z-scored BOLD observations. Since our goal is to derive the (dynamic) edge-centric properties from the (static) nFCterm·mth matrix ( R = [ r i j ]), we need to define a null hypothesis that discounts any temporal dependencies but retains the observed spatial correlations in R . A simple null hypothesis on the distribution of (the columns of) z that satisfies this criterion is Z ( t ) ~ N ( 0 , R ) , that is, i.i.d. multivariate Gaussian RV with a covariance matrix R matching the observed nFCterm·mth. If we denote the state of the system at time t as the column vector z ( t ) = [ z 1 ( t ) , … , z N ( t ) ] ⊤ , the null hypothesis simply states that z ( t ) is drawn from the same multivariate Gaussian distribution at each time t , independently of the other samples.

Derivation of edge FCterm·mth With capital letters denoting RVs, the expected product between two edge time seriesterm·key C j k ( t ) and C l m ( t ) at time t is 8 E [ C j k ( t ) C l m ( t ) ] = E [ Z j ( t ) Z k ( t ) Z l ( t ) Z m ( t ) ] = κ ( Z j ( t ) Z k ( t ) Z l ( t ) Z m ( t ) ) + E [ Z j ( t ) Z k ( t ) ] E [ Z l ( t ) Z m ( t ) ] + E [ Z j ( t ) Z l ( t ) ] E [ Z k ( t ) Z m ( t ) ] + E [ Z j ( t ) Z m ( t ) ] E [ Z k ( t ) Z l ( t ) ] , which follows from the definition of the joint cumulant κ ( Z j ( t ) Z k ( t ) Z l ( t ) Z m ( t )), also noting that the products involving the expectation of a single variable are equal to zero (i.e., E [ Z i ( t ) ] = 0 ) since z i are z-scored. The expression of the expectation in terms of joint cumulants is sometimes referred to as moment-cumulants formula 38 , 39 . The joint cumulant is equal to zero for Gaussian RVs 39 , allowing a simplification of Eq. ( 8 ) known as Isserlis’ theorem 40 : 9 E [ C j k ( t ) C l m ( t ) ] = E [ Z j ( t ) Z k ( t ) ] E [ Z l ( t ) Z m ( t ) ] + E [ Z j ( t ) Z l ( t ) ] E [ Z k ( t ) Z m ( t ) ] + E [ Z j ( t ) Z m ( t ) ] E [ Z k ( t ) Z l ( t ) ] . If we additionally assume that Z ( t ) is ergodic (which does not preclude it from being a multivariate autoregressive process 19 ), all the involved terms become independent of time and Eq. ( 9 ) further simplifies to 10 E [ C j k ( t ) C l m ( t ) ] = r j k r l m + r j l r k m + r j m r k l . In a mathematical sense, the expectation is to be intended over the population (ensemble), whereas the eFCterm·mth is computed from a single session (sample). However, the ergodic assumption guarantees that sample estimates converge to the ensemble expectation as the number of time frames increases. We can then obtain an estimate of the eFCterm·mth by substituting Eq. ( 10 ) into Eq. ( 7 ): 11 eFCterm·mth j k , l m = r j k r l m + r j l r k m + r j m r k l 1 + 2 r j k 2 1 + 2 r l m 2 . The i.i.d. multivariate Gaussian null modelterm·key clearly satisfies the ergodic hypothesis since it has no memory. However, note that neither the i.i.d. nor the Gaussian properties are necessary since the derivations in this section assume ergodicity with the only additional constraint that the joint cumulant is equal to zero.

Derivation of edge communities In order to formalise the intuition that two edges ( j k ) and ( j ′ k ′ ) with similar rows of the eFCterm·mth are likely to be clustered together, let us define their distance d j k , j ′ k ′ as the ℓ 1 norm of the difference between the corresponding rows of the (unnormalised) eFCterm·mth: 12 d j k , j ′ k ′ : = ∑ l , m = 1 N 1 T ∣ E C j k C l m ⊤ − E C j ′ k ′ C l m ⊤ ∣ = ∑ l , m = 1 N ∣ r j k r l m + r j l r k m + r j m r k l − r j ′ k ′ r l m − r j ′ l r k ′ m − r j ′ m r k ′ l ∣ = ∑ l , m = 1 N ∣ z j z k ⊤ z l z m ⊤ + z j z l ⊤ z k z m ⊤ + z j z m ⊤ z k z l ⊤ − z j ′ z k ′ ⊤ z l z m ⊤ − z j ′ z l ⊤ z k ′ z m ⊤ − z j ′ z m ⊤ z k ′ z l ⊤ ∣ = ∑ l , m = 1 N 3 ∣ z j ( z l ⊤ z m ) z k ⊤ − z j ′ ( z l ⊤ z m ) z k ′ ⊤ ∣ . We now have the necessary ingredients to build a measure of similarity between the nodes, which can be used to predict the edge cluster similarity in ref. 17 . There, the similarity between two nodes is measured as the frequency with which the corresponding edges are clustered together (having fixed the number of communities to 10). Instead of discrete assignments to 10 communities, Eq. ( 12 ) provides a continuous measure of the distance between two edges. The distance between two nodes i and j can then be defined as the sum of the distances between the edges starting from i and j and reaching the same target : 13 d i , j = ∑ k = 1 N d i k , j k = ∑ k , l , m 3 ∣ ( z i − z j ) ( z l ⊤ z m ) z k ⊤ ∣ ≤ ( z i − z j ) ∑ k , l , m 3 ( z l ⊤ z m ) z k ⊤ = ∑ t ( z i ( t ) − z j ( t ) ) 2 1 2 ∑ k , l , m 3 ( z l ⊤ z m ) z k ⊤ = ( T − 1 ) ( Var [ z i ] + Var [ z j ] − 2 r i j ) 1 2 ∑ k , l , m 3 ( z l ⊤ z m ) z k ⊤ = 1 − r i j 1 2 ( 2 ( T − 1 ) ) 1 2 ∑ k , l , m 3 ( z l ⊤ z m ) z k ⊤ ∝ 1 − r i j 1 2 , where we used the Cauchy-Schwarz inequality and noted that the terms in the square brackets form a constant (independent of i and j ). It is then apparent that the edge-cluster similarity 17 between nodes i and j can be approximated by the nFCterm·mth. Once again, note that the i.i.d. Gaussian RV assumption can be relaxed since the derivations in this section are based on Eq. ( 10 ), which requires ergodicity with the only additional constraint that the joint cumulant is equal to zero.

Derivation of RSSterm·mth from the BOLD signal Recalling the definition of the edge time seriesterm·key c i j ( t ) in Eq. ( 6 ), the RSSterm·mth defined in ref. 16 can be approximated as the squared Euclidean norm of the z-scored BOLD signal, up to a constant factor: 14 RSSterm·mth ( t ) : = ∑ i < j c i j ( t ) 2 = 1 2 ∑ i , j = 1 N c i j ( t ) 2 − ∑ i = 1 N c i i ( t ) 2 = 1 2 ∑ i , j = 1 N z i ( t ) 2 z j ( t ) 2 − ∑ i = 1 N z i ( t ) 4 = 1 2 z ( t ) 4 − ∑ i = 1 N z i ( t ) 4 ≈ 1 2 z ( t ) 2 . The approximation does not rely on the i.i.d. assumption; it is valid under the Gaussian null hypothesis and, more generally, for distributions with finite kurtosis—including fMRI data. Under this assumption, z ( t ) 4 dominates ∑ i z i ( t ) 4 in Eq. ( 14 ), as can be seen from the ratio of their (expected) values: 15 E ∑ i Z i ( t ) 4 E Z ( t ) 4 ≤ ∑ i Kurt [ Z i ( t ) ] N 2 → N → ∞ 0 .

The approximation in Eq. ( 14 ) can be replaced by an exact equality if all the N 2 edge time seriesterm·key are included in the RSSterm·mth definition (that is, all the ( i , j ) tuples, rather than only the pairs with i < j ): 16 RSSterm·mth all ( t ) : = ∑ i , j c i j ( t ) 2 = z ( t ) 2 .

Why are frames with the largest RSSterm·mth most similar to the nFCterm·mth?

Having rewritten the RSSterm·mth as the squared Euclidean norm in Eq. ( 14 ), we can more easily investigate the conditions underpinning the largest RSSterm·mth fluctuations. Let us withen the BOLD vector Z ( t ) to obtain the RV 17 W ( t ) : = R − 1 2 Z ( t ) and let 18 R = U Λ U ⊤ be the eigendecomposition of the correlation matrix R , with Λ = diag( λ 1 , …, λ N ) and U being the unitary matrix of eigenvectors of R (since R is symmetric). Without loss of generality, assume that the eigenvaluesterm·key are sorted in descending order, such that λ 1 is the largest eigenvalueterm·key and u 1 is the corresponding leading eigenvector. The RSSterm·mth can be treated as a RV and rewritten in terms of W ( t ) and the eigenvector matrix U : 19 RSSterm·mth ( t ) ≈ 1 2 Z ( t ) 2 = 1 2 Z ( t ) ⊤ Z ( t ) = 1 2 W ( t ) ⊤ ( R 1 2 ) ⊤ R 1 2 W ( t ) = 1 2 ( W ( t ) ⊤ U ) Λ ( U ⊤ W ( t ) ) = 1 2 ∑ i λ i ( U ⊤ W ( t ) ) i 2 = 1 2 ∑ i λ i u i , W ( t ) 2 20 = 1 2 ∑ i λ i u i 2 W ( t ) 2 cos 2 Θ i ( t ) , where u i is the i -th eigenvector and Θ i ( t ) is the RV representing the angle formed by the vectors u i and W ( t ) at time t . Also note that u i 2 = 1 because U is unitary. For any realisations w ( t ) with squared norm w ( t ) 2 , an upper bound on the RSSterm·mth is obtained as 21 RSSterm·mth ( t ) ≤ 1 2 w ( t ) 2 max i λ i ∑ i cos 2 θ i ( t ) = 1 2 λ max w ( t ) 2 , noting that 22 ∑ i cos 2 θ i ( t ) = ∑ i ⟨ u i , w ( t ) ⟩ 2 u i 2 w ( t ) 2 = U w ( t ) 2 w ( t ) 2 = 1 . The upper bound is reached when θ 1 ( t ′ ) = 0 or θ 1 ( t ′ ) = π , which implies that w ( t ′ ) = c u 1 , where c is a constant. In other words, w ( t ′ ) is aligned (parallel or antiparallel) with the leading eigenvector u 1 . When this happens, the BOLD signal vector z ( t ′ ) must also be aligned with u 1 : 23 z ( t ′ ) = R 1 2 w ( t ′ ) = R 1 2 c u 1 = U Λ 1 2 U ⊤ u 1 = c λ 1 1 2 u 1 . We can then refine our theoretical understanding of the RSSterm·mth peaks: not only they occur when the Euclidean norm of the BOLD signal vector is large (as per Eq. ( 14 )) but, most likely, when it is well aligned with the leading eigenvector of the static nFCterm·mth (see Fig. 5 b). If the alignment were perfect at a frame t ′ , the instantaneous estimate of the nFCterm·mth would be 24 z ( t ′ ) z ( t ′ ) ⊤ = c 2 λ 1 u 1 u 1 ⊤ , i.e., an approximation of the nFCterm·mth obtained from its leading eigenvector only (and independent of the sign of c ). This approximation would achieve an average similarity (Pearson correlation coefficient) of r = 0.69 with the nFCterm·mth over 100 unrelated participants of the HCP dataset (while, in practice, the highest similarity achieved by the top frame was r = 0.53). However, the alignment with u 1 need not be perfect: in general, large RSSterm·mth values can be expected whenever the BOLD signal is a mixture of the top eigenvectors. Additional principal components are expressed as more large-RSS frames are averaged, suggesting why the top 5% frames alone are sufficient for an almost perfect reconstruction of the nFCterm·mth (Fig. 3 ). We have thus explained why cofluctuationterm·key patterns corresponding to frames with the largest RSSterm·mth exhibit the highest similarity with the nFCterm·mth. These results are based on Eq. ( 14 ) and hold true under the assumption of finite kurtosis (which also applies in the specific case of the null hypothesis, i.e., for Gaussian variables). The i.i.d. assumption is not required.

Null distribution of the RSSterm·mth The RSSterm·mth can be written as a simple quadratic form RSSterm·mth ( t ) = 1 2 Z ( t ) 2 = 1 2 Z ( t ) ⊤ Z ( t ) , which is known to follow a generalised χ 2 distribution under the null hypothesis of Gaussian variables 41 . The weights of the non-central chi-square components are proportional to the eigenvaluesterm·key of the nFCterm·mth matrix, i.e., λ 1 2 , … , λ N 2 . Another characterisation of this distribution is provided by Eq. ( 19 ): under the null hypothesis, the inner product 〈 u i , W ( t )〉 follows a normal Gaussian distribution since W ( t ) ~ N ( 0 , 1 ) and U is unitary. Therefore, ⟨ u i , W ( t ) ⟩ 2 follows a χ 2 distribution and each term λ i 2 ⟨ u i , W ( t ) ⟩ 2 in Eq. ( 19 ) follows a Gamma ( k = 1 2 , θ = 2 λ i ) distribution. The RSSterm·mth is thus obtained as a sum of N independent Gamma-distributed RVs, each associated with one eigenvalueterm·key of the nFCterm·mth. The tail of the RSSterm·mth is best approximated by the RVs associated with the largest eigenvaluesterm·key (which have the largest mean and variance), while including smaller eigenvaluesterm·key provides an increasingly fuller characterisation of the whole distribution (Fig. 4 b). The mean and variance of the RSSterm·mth can be readily obtained from the properties of the Gamma distribution: 25 E [ RSSterm·mth ] = 1 2 ∑ i λ i = N 2 26 Var [ RSSterm·mth ] = ∑ i λ i 2 . Higher moments of the RSSterm·mth null distribution can be derived from its moment-generating function: 27 M RSSterm·mth ( s ) = ∏ i ( 1 − 2 λ i s ) − 1 2 .

Why are large RSSterm·mth fluctuations present in many datasets?

The moment-generating function in Eq. ( 27 ) can be employed to show that the RSSterm·mth is subexponential under the null hypothesis, which explains its heavy tail and the consequent large events 25 , 42 . Specifically, the subexponential feature of the null RSSterm·mth follows from the sufficient condition 28 M RSSterm·mth − E [ RSSterm·mth ] ( s ) = ∏ i ( 1 − 2 λ i s ) − 1 2 exp − λ i s 2 ≤ ∏ i exp λ i 2 s 2 = exp s 2 ∑ i λ i 2 , ∀ ∣ s ∣ ≤ ( 4 λ max ) − 1 . However, we can expect this behaviour under the more general hypothesis that the z-scored BOLD signal is sub-Gaussian, i.e., its tail decays at least as fast as that of a Gaussian RV (including, for example, any uniformly-bounded RVs). The reason is that the square of a sub-Gaussian RV is subexponential, and the sum of independent subexponential RVs is also subexponential. Therefore, being the RSSterm·mth closely approximated by a sum of squared RVs (as per Eq. ( 14 )), extreme events are to be expected under the general sub-Gaussian assumption for the BOLD signal, which offers an explanation for the large RSSterm·mth fluctuations observed in most fMRI datasets.

How do nFCterm·mth modules influence the edge cofluctuationsterm·key?

Interestingly, Pope et al. 26 have recently reported a connection between the presence of structural modules and the occurrence of large events in the edge cofluctuationsterm·key (RSSterm·mth). Insofar as structural and functional modules are in agreement 43 , 44 , we can explain these findings based on the nFCterm·mth spectrum. How do functional modules shape the eigenspectrum of the nFCterm·mth? In the ideal case of a block-diagonal matrix (with zeroes outside the blocks), the sum of the eigenvaluesterm·key corresponding to each block coincides with the block size (since the diagonal elements are all ones and the trace is preserved under diagonalisation). As such, the largest eigenvalueterm·key is bounded by the size of the largest block, i.e., larger functional modules allow for larger eigenvaluesterm·key. In turn, large eigenvaluesterm·key underpin the high-amplitude cofluctuationsterm·dom, as shown in the ‘Why are frames with the largest RSSterm·mth most similar to the nFCterm·mth?’ section. Therefore, if the size of the modules is reduced via randomisation of the structural connectivity as in ref. 26 (SI Fig. 3) , the expected magnitude of the RSSterm·mth cofluctuationsterm·key will drop according to Eq. ( 21 ). This offers a mathematical explanation for the lower RSSterm·mth event count when the modular structure is disrupted.

A null modelterm·key for binary edge time seriesterm·key Consider two parcels j and k with z-scored BOLD activity Z j ( t ) and Z k ( t ) and correlation coefficient r j k as defined by the nFCterm·mth matrix. Under the static Gaussian null hypothesis, Z k ( t ) ~ N ( 0 , 1 ) and 29 Z j ( t ) = r k j Z k ( t ) + 1 − r k j 2 V j , where V j ~ N ( 0 , 1 ) , and the square root coefficient ensures unit variance for Z j ( t ) since both parcels are z-scored. The edge time seriesterm·key corresponding to the two parcels is C j k ( t ) = Z j ( t ) Z k ( t ), as per Eq. ( 6 ). By definition, the binary edge time seriesterm·key C j k ( t ) ¯ is equal to 1 when the cofluctuationsterm·key between j and k are positive, i.e., when Z j ( t ) and Z k ( t ) have the same sign. In other words, C j k ( t ) ¯ is a Bernoulli( p j k ) RV with probability 30 p j k = P [ C j k ( t ) ¯ > 0 ] = P [ Z j ( t ) Z k ( t ) > 0 ] = 2 P [ ( Z k ( t ) > 0 ) ∩ ( Z j ( t ) > 0 ) ] = 2 P ( Z k ( t ) > 0 ) ∩ ( r k j Z k ( t ) + 1 − r k j 2 V j > 0 ) = 2 P ( Z k ( t ) > 0 ) ∩ V j Z k ( t ) > − r k j 1 − r k j 2 = 2 P [ A ] , where A is the event ( Z k ( t ) > 0 ) ∩ V j Z k ( t ) > − r k j 1 − r k j 2 . From a geometric perspective, A is satisfied by any vector ( Z k , V j ) whose polar angle is between π /2 and arctan ( − r j k 1 − r k j 2 ) . Therefore, 31 p j k = 2 P [ A ] = 2 π 2 − arctan ( − r j k 1 − r k j 2 ) 2 π = 1 2 + arcsin ( r j k ) π . In conclusion, under the null hypothesis, the binary edge time seriesterm·key C j k ( t ) ¯ can be modelled as a Bernoulli( p j k ) RV with p j k = 1 2 + arcsin ( r j k ) π .

Relationship with coactivation patternsterm·mth (CAPsterm·mth) Here, we will focus on explaining a well-known finding in the CAPsterm·mth literature: choosing a seed node, Liu and Duyn 12 observed a strong correlation between the BOLD activity at frames where the seed is highly active and the corresponding seed-based FCterm·mth. Note that the correlation is measured between two N -dimensional vectors (where N is the number of nodes): the first vector is the BOLD signal and the second vector is the column of the nFCterm·mth matrix that corresponds to the chosen seed node. In the following, we will show that the i.i.d. Gaussian null modelterm·key can predict the observed correlation. Let k be the seed node and z k its z-scored BOLD time series. As usual, the associated random process is denoted with the capital letter Z k ( t ). If we condition on a specific value of the seed node, say z k * : = Z k ( t * ) , the BOLD activity vector at the corresponding time t * is expected to align with the k -th column of the nFCterm·mth (denoted as R ^ k ): 32 E [ Z ( t * ) ] = E [ Z ( t ) ∣ Z k ( t ) = z k * ] = z k * R ^ k . This is an elementary property of multivariate Gaussian RVs that follows directly from Eq. ( 29 ): 33 E [ Z j ( t * ) ] = E [ Z j ( t ) ∣ Z k ( t ) = z k * ] = E r k j Z k ( t ) + 1 − r k j 2 V j ∣ Z k ( t ) = z k * = z k * r k j , ∀ j = 1 , … , N . Perhaps less intuitively, the conditional correlation between the BOLD vector Z ( t * ) and the k -th column of the sample nFCterm·mth ( R ^ k ) is also proportional to the seed value z k * . In order to prove it, let us first compute their expected inner product: 34 E [ ⟨ Z ( t * ) , R ^ k ⟩ ] = E ∑ j N Z j ( t * ) 1 T − 1 ∑ t T Z j ( t ) Z k ( t ) ≈ ∑ j E Z j ( t * ) 1 T − 1 ∑ t ≠ t * E Z j ( t ) Z k ( t ) = ∑ j z k * r k j 1 T − 1 ∑ t ≠ t * r k j = z k * R k 2 . Similarly, the covariance can be approximated as 35 Cov [ Z ( t * ) , R ^ k ] = E ⟨ Z ( t * ) , R ^ k ⟩ N − 1 − ∑ j Z j ( t * ) ∑ l r ^ k l N ( N − 1 ) ≈ z k * R k 2 N − ∑ j ≠ l r k j r k l N ( N − 1 ) The expected sample covariance is directly proportional to z k * and thus peaks at frames with the highest activity of the seed node k . This remains true after normalising the covariance to obtain the correlation coefficient, as shown in Fig. 6 a.

Human Connectome Projectterm·key fMRI dataset This study used openly-available and independently-acquired resting-state fMRI (rsfMRI) data from the Human Connectome Projectterm·key (HCP) S1200 release 45 . In particular, we used the “100 unrelated subjects” dataset: a subset of 100 non-twins adult participants which were pre-selected by the HCP coordinators (54% female; mean age = 29.11 ± 3.67 years; age range, 22–36 years). The HCP study was approved by the Washington University Institutional Review Board, and informed consent was obtained from all participants. All subjects were scanned on a customized Siemens 3T “Connectome Skyra” with a 32-channel head coil, housed at Washington University in St. Louis. rsfMRI data were acquired in four runs of 15 min over a 2-day period, with eyes open and relaxed fixation on a projected bright cross-hair on a dark background (presented in a darkened room). Resting-state images were collected with the following parameters: gradient-echo EPI sequence, run duration = 14:33 min, TR = 720 ms, TE = 33.1 ms, flip angle = 52 ∘ , FOV = 208 × 180 mm (RO x PE), matrix = 104 × 90 (RO x PE), slice thickness = 2 mm, 2-mm isotropic voxel resolution, multi-band factor = 8, echo spacing = 0.58 ms, BW = 2290 Hz/Px).

Preprocessing and ICA-FIX denoising Functional images in the HCP dataset were minimally preprocessed according to the pipeline described in ref. 46 . In short, the data were corrected for gradient distortion, susceptibility distortion and motion and then aligned to a corresponding T1-weighted image with one spline interpolation step. This volume was further corrected for intensity bias, normalised to a mean of 10,000, projected to the 32k_fs_LR mesh (excluding outliers), and aligned to a common space using a multi-modal surface registration.

In addition, the preprocessed rsfMRI data were cleaned of structured noise through a process that pairs independent component analysis (MELODIC) with FIX to automatically remove non-neural spatiotemporal components (trained on 25 hand-labelled HCP subjects). The FIX approach and initial results of classification accuracy are detailed in ref. 47 , and the effects of the ICA + FIX cleanup (and optimal methods to remove the artefactual components from the data) are evaluated in detail in ref. 48 . The cleaning pipeline is described more comprehensively in the HCP S1200 release reference manual ( https://humanconnectome.org/study/hcp-young-adult/document/1200-subjects-data-release/ ) and the preprocessing and the cleaning scripts are openly available on Github ( https://github.com/Washington-University/HCPpipelines/ ). The resulting ICA-FIX denoised rfMRI grayordinate surface time series are available as CIFTI files following the naming pattern: *REST1,2_LR,RL_Atlas_MSMAll_hp2000_clean.dtseries.nii.

The Schaefer200 parcellation was used to define 200 areas on the cerebral cortex 49 . This functional parcellation was designed to optimise both local gradient and global similarity measures of the fMRI signal and is openly available in ‘32k fs LR’ space for the HCP dataset. The nodes are mapped to the Yeo canonical functional networks 50 . The parcellated data were analysed both before and after regressing the global signal. The theoretical derivations and predictions hold and perform equally well in both cases, and we report any significant differences when they occur. Unless otherwise stated, the GSR results are shown in the figures since they are more directly comparable to those published in refs. 16 , 17 , noting in particular that GSR was performed in ref. 16 . Despite the ICA-FIX preprocessing pipeline used here is entirely different from those employed in refs. 16 , 17 , our results are in excellent agreement with the previously published ones.

Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Supplementary information Supplementary Information Reporting Summary