Abstract
Individual differences in the spatial organization of resting-state networks have received increased attention in recent years. Measures of individual-specific spatial organization of brain networks and overlapping network organization have been linked to important behavioral and clinical traits and are therefore potential biomarker targets for personalized psychiatry approaches. To better understand individual-specific spatial brain organization, this paper addressed three key goals. First, we determined whether it is possible to reliably estimate weighted (non-binarized) resting-state network maps using data from only a single individual, while also maintaining maximum spatial correspondence across individuals. Second, we determined the degree of spatial overlap between distinct networks, using test-retest and twin data. Third, we systematically tested multiple hypotheses (spatial mixing, temporal switching, and coupling) as candidate explanations for why networks overlap spatially. To estimate weighted network organization, we adopt the Probabilistic Functional Modes (PROFUMO) algorithm, which implements a Bayesian framework with hemodynamic and connectivity priors to supplement optimization for spatial sparsity/independence. Our findings showed that replicable individual-specific estimates of weighted resting-state networks can be derived using high-quality fMRI data within individual subjects. Network organization estimates using only data from each individual subject closely resembled group-informed network estimates (which was not explicitly modeled in our individual-specific analyses), suggesting that cross-subject correspondence was largely maintained. Furthermore, our results confirmed the presence of spatial overlap in network organization, which was replicable across sessions within individuals and in monozygotic twin pairs. Intriguingly, our findings provide evidence that overlap between 2-network pairs is indicative of coupling. These results suggest that regions of network overlap concurrently process information from both contributing networks, potentially pointing to the role of overlapping network organization in the integration of information across multiple brain systems.
1 Introduction
Recent studies have revealed substantial inter-individual variability in the spatial organization of the brain as measured with resting state functional MRI (rfMRI) (Braga & Buckner, 2017; D. Wang et al., 2015; Glasser et al., 2016; Gordon, Laumann, Adeyemo, & Petersen, 2017; Gordon, Laumann, Adeyemo, Gilmore, et al., 2017; Gordon, Laumann, Gilmore, et al., 2017; Laumann et al., 2015). Importantly, such inter-individual spatial variability in functional brain organization is strongly associated with behavioral traits (Bijsterbosch et al., 2018; Kong et al., 2019). The overarching objective of this paper is to characterize weighted (i.e., non-binarized) spatial organization of resting-state networks within individuals, with a specific focus on gaining insights into spatially overlapping network organization. We focus on network overlap because recent work has highlighted that individual differences in network overlap are strongly associated with behavior (Bijsterbosch et al., 2019). Yet, the complex spatially overlapping nature of network organization is largely underestimated and understudied in traditional rfMRI analysis approaches that enforce binary area/network assignment or encourage spatial independence of networks.
Identifying network organization at the individual level rather than using group information raises multiple challenges. First, individual estimates of network organization are noisier than group estimates, partly because group estimates benefit from collating large amounts of data. This challenge can be partially addressed by obtaining large amounts of data from each individual (precision functional mapping approach (Gratton et al., 2020)), but such extensive data acquisition may not be feasible in all participants and settings. Second, although a purely individual-specific set of network maps represents the most accurate and unbiased estimate of the individual’s brain organization (Gordon, Laumann, Gilmore, et al., 2017), it may lack network correspondence across individuals. Assuming the presence of cross-participant commonalities in their network structure, such group correspondence is valuable for network labeling and interpretability, and essential for group-level and between-subject analytical comparisons. Group-based estimates applied to individuals have built-in correspondence, but these individual estimates may be biased towards the group estimate (Bijsterbosch et al., 2019). Probabilistic Functional Modes (PROFUMO) is a Hierarchical Bayesian algorithm developed to try to optimize this trade-off by using group-level priors to achieve correspondence, whilst optimizing individual-specific estimates to maximize the accuracy of individual network maps (Harrison et al., 2015). Compared to other hierarchical network approaches such as Kong et al. (2019), PROFUMO estimates the group prior from the data (instead of adopting an atlas) and includes a group prior on the temporal connectivity matrix (for a comprehensive overview of hierarchical model differences, please see Appendix 1 for Farahibozorg et al., 2021). Although PROFUMO has been successfully applied in group data (Bijsterbosch et al., 2018; Farahibozorg et al., 2021; Harrison et al., 2020), an open question is whether it can be robustly applied to data from only a single individual without sacrificing correspondence, which is of particular interest in the context of personalized psychiatry and translational work in non-human primates and other animal models. The first goal of this paper was to determine whether PROFUMO can reliably estimate weighted (i.e., non-binarized) resting-state networks using only data from a single subject, while also achieving substantial “non-enforced” correspondence across individuals (i.e., without data from other participants contributing to a group prior to inform correspondence).
Spatial overlap between rfMRI networks beyond classical “hub regions” has been observed across a variety of analytical brain representations (Lee et al., 2016; Lin et al., 2018; Najafi et al., 2016), and has been linked to behavioral traits (Bijsterbosch et al., 2019). PROFUMO accurately estimates spatial overlap in rfMRI network organization (Bijsterbosch et al., 2019), which is a key advantage compared to approaches that aim for a “hard” binarized parcellation or approaches that enforce spatial independence between networks (Bijsterbosch et al., 2020). Despite broad interest in “hub” regions within the graph theory domain that typically adopts hard parcellations (Bertolero et al., 2018; Buckner et al., 2009; Warren et al., 2014), a detailed spatial investigation into individual-specific overlap between weighted estimates of network organization is lacking. Studying the overlapping properties of individualized brain networks is of interest because network overlap may be a sensitive marker for use in personalized psychiatry settings (Gratton et al., 2020; Insel, 2014; Williams, 2016) given prior evidence of behavioral relevance (Bijsterbosch et al., 2019), provided that it can be robustly and reliably detected in individuals. Furthermore, this work contributes to the broader literature on precision functional mapping (Gordon, Laumann, Gilmore, et al., 2017; Gratton et al., 2018; Poldrack, 2017; Poldrack et al., 2015). The second goal of this paper was to characterize network overlap among weighted individual-specific resting state networks estimated using PROFUMO.
Each type of resting-state fMRI analysis (e.g., using different approaches such as Independent Component Analysis, graph theory, amplitude of low frequency fluctuations, PROFUMO, etc.) provides a different low-dimensional representation of the dataset (Bijsterbosch et al., 2020). Although optimized to best fit the data (Bijsterbosch et al., 2021), these brain representations (i.e., summary measures derived from different approaches to rfMRI data, such as connectivity matrices, network maps, amplitude maps, etc.) are necessarily lossy given the intrinsic goal of dimensionality reduction. In the case of PROFUMO, a set of stationary large-scale modes of brain organization are derived that collapse fine-grained spatial structure and simplify temporally dynamic processes. As such, there are multiple potential mechanisms that may give rise to spatially overlapping network organization as observed between network maps estimated with PROFUMO. First, it is possible that a brain region in which two networks appear to overlap may, in fact, be a spatially heterogeneous mixture of cortical patches that are individually part of either network 1 or network 2, implying no real functional “link” between the two networks as a result of the overlap (Fig. 1A). Such a spatially heterogeneous overlap region may, for example, result from network interdigitation (Braga & Buckner, 2017; Braga et al., 2019), or regional gradients (Blazquez Freches et al., 2020; Haak et al., 2018). Second, the region of network overlap may dynamically switch its network allegiance over time to be part of either network 1 or network 2 at any given moment (Hutchison et al., 2013) (Fig. 1B). Such dynamic switching would appear as a spatially overlapping network structure given the stationary (time-averaging) nature of the PROFUMO model. Third, network overlap might indicate that signals from network 1 and network 2 are jointly processed and “deeply functionally integrated” within regions of network overlap (Fig. 1C). This third hypothesis is perhaps the most intriguing as it may indicate information coupling involving a specific functional role of overlap regions, and contributing to between-network communication. The third goal of this paper was to systematically compare these spatial mixture, dynamic switching, and coupling hypotheses of network overlap.
In this paper, we leveraged a unique subset of the Human Connectome Project Young Adult data (Van Essen et al., 2013), by focusing on individuals who underwent three complete resting-state sessions (3 T, 3 T re-scan, and 7 T) for a total of approximately 3 hours of rfMRI data per individual. The resulting sample of N = 20 that met this criteria further included 8 monozygotic twin pairs, providing a rich cohort to investigate individual-specific weighted network organization. Our focus on a small sample of densely sampled individuals was informed by the interest in individual-specific network organization. Although no brain-behavior associations were feasible given the small sample (Marek et al., 2022), prior work has extensively studied individual differences in PROFUMO, including behavioral associations with spatial organization (Bijsterbosch et al., 2018), heritability (Harrison et al., 2020), and network variants as a function of dimensionality (Farahibozorg et al., 2021). The results of this work support the application of PROFUMO weighted networks in individual participants, which paves the way for future applications in a personalized psychiatry framework. Furthermore, our findings suggest a coupling mechanism underlying spatially overlapping network organization, which supports the hypothesis that regions of network overlap play an important functional role in terms of cross-network coupling (Gordon et al., 2018).
2 Methods
2.1 Dataset
We used high-quality data from the Human Connectome Project (HCP) (Van Essen et al., 2013), focusing on N = 20 individuals (including 8 monozygotic twin pairs) who underwent a complete set of four 3 T, four 7 T, and four retest 3 T runs, thereby accumulating approximately 3 hours (13,200 timepoints) of rfMRI data across 12 scans per individual. This HCP-YA sub-sample was 80% female with a mean age of 30.1 years (standard deviation = 3.84; range 22-34). Briefly, the 3 T rfMRI data were acquired at 2 mm isotropic voxel size using a multiband factor of 8, a TR of 0.72 seconds, and a TE of 33 ms; the 7 T rfMRI were acquired at 1.6 mm isotropic voxel size, a multiband acceleration of 5, in-plane acceleration 2, a TR of 1.0 seconds, and a TE of 22.2 ms (see further details in Smith et al., 2013; Van Essen et al., 2013; Vu et al., 2017). Data were preprocessed using the HCP minimal processing pipelines (Glasser et al., 2013). All 3 T and 7 T data were analyzed in the Connectivity Informatics Technology Initiative (CIFTI) format, which consists of 91,282 grayordinates with approximately 2-mm spatial resolution (i.e., 7 T data were downsampled to match the 3 T spatial resolution). ICA-FIX was then applied to remove structured noise (Griffanti et al., 2014; Salimi-Khorshidi et al., 2014), and data were aligned using multimodal surface matching (MSM-All; Robinson et al., 2014) to align areal features (myelin and RSNs).
2.2 PROFUMO estimation
PROFUMO is a matrix factorization approach for the estimation of resting-state networks that adopts spatial priors, temporal priors, and a noise model in a hierarchical Bayesian model (see Harrison et al., 2015, 2020 and Supplementary Tables S1 and S2). PROFUMO spatial maps are calculated by vertex-wise multiplying the probability (0-1; the probability that a given weight is drawn from the signal rather than the noise distribution) by the estimated mean (derived from Gaussian mixture model). PROFUMO was applied in several distinct ways:
Classic group-PROFUMO was performed in which data from all 12 runs across all 20 participants were used and modeled hierarchically according to the levels of subjects and (beneath that) runs. Notably, this version of PROFUMO is recommended for studies that include only one or a small number of sessions per individual.
Single-subject PROFUMO was performed independently for each of the 20 participants, using all 12 runs for each, those 12 being considered as separate “subjects” in the estimation of PROFUMO’s Gaussian mixture model (Harrison et al., 2020). As such, the “group” priors for single-subject PROFUMO here (and in #4 & #5 below) reflect the single-subject averages, such that no data/information from other participants was used in these analyses.
Classic group-PROFUMO was performed similar to case (#1), using 12 individual runs from 12 separate participants. This was in order to obtain a group-level estimation of modes using the same amount of data available for the separate subject runs (case (#2), i.e., matching the effective signal to noise ratio; SNR).
Test-retest single-subject PROFUMO was performed independently for each of the 20 participants and independently using two sets of 6 runs each (split evenly across 3 T, 7 T, and retest data; and always including a pair of opposite phase-encode directions).
Single-subject PROFUMO similar to case (#2) was performed stepwise on cut-down versions of the 12 runs that included a progressively increasing number of timepoints (in increments of 1/12th of the run timepoints) to determine how much data are needed to obtain reliable network estimates. This approach was taken to ensure that all PROFUMO runs included equal contributions from 3 T, retest, and 7 T data across all phase-encode directions.
All of the PROFUMO runs were performed at a dimensionality of 20 to focus on the spatial organization of large-scale resting-state networks. Furthermore, classic group-PROFUMO (#1) was repeated at dimensionalities 15, 30, 40, and 50 to determine the stability of the selected networks across PROFUMO decomposition dimensionalities. Beyond the key parameter of dimensionality (number of networks), all PROFUMO parameters were set to the default. A summary of all PROFUMO parameters, hyperparameters, and hyperpriors can be seen in Supplementary Tables S1 and S2. Beyond these PROFUMO parameters, additional key measures of interest for this paper include the similarity between group and individual network estimates indicative of correspondence and the test-retest reliability of network spatial maps representative of stable and reliable network estimates (described further below in Sections 2.3 and 2.4).
2.3 PROFUMO mode selection
The Hungarian algorithm (a.k.a. “munkres” algorithm) was used to reorder PROFUMO modes (i.e., “networks”) for each of the single-subject runs (#2 above), the split 1 and split 2 single-subject runs (#4 above), and the 12-run group run (#3 above) to best-match the mode order obtained from the full group results. Briefly, the Hungarian algorithm solves the assignment problem by permuting rows of the network-to-network spatial correlation matrix to minimize the trace of the permuted cost matrix (Kuhn, 1955; Munkres, 1957). For each of the 20 modes, we estimated the test-retest correlation per subject as the Pearson correlation across all 91,282 CIFTI grayordinates between split 1 and split 2 single-subject runs. For each of the 20 modes, we also estimated the subject-group correlation per subject as the Pearson correlation across all 91,282 vertices between the map from the subject run (#2 above) and the subject-specific estimated map from the group run (#1 above). This “group-individual run” measure addresses the question of correspondence across individuals because it compares classic group PROFUMO modes (which benefit from explicit correspondence through the group prior in the hierarchical Bayesian algorithm) to single-subject PROFUMO modes (which do not involve any group information). Modes that achieved both a median (across subjects) test-retest correlation of 0.7 or greater and a median (across subjects) group-individual correlation of 0.7 or greater were used in subsequent analyses. Notably, 0.7 was previously recommended as a default lowest acceptable standard of reliability in basic research (Nunnally, 1967). Out of the 20 PROFUMO modes, 12 modes met these requirements. Within these 12 modes, individual participants’ missing modes were defined as modes with both a test-retest correlation lower than 0.2 and a subject-group correlation lower than 0.2 (both estimated at the subject level). Missing modes were ignored in subsequent analyses. For naming purposes, modes were spatially mapped onto the Yeo-7 parcellation (Yeo et al., 2011) and we followed the network naming taxonomy suggested by Uddin et al. (2019).
2.4 PROFUMO mode stability
We assessed the stability of spatial maps in several ways. As described in Section 2.3, the first two measures of mode stability were within-participant test-retest reliability between splits 1 and 2, and within-subject similarity (group-individual run) correlations between the subject run (#2 above) and the subject estimates obtained in the group run (#1 above). We also tested twin similarity (individual runs) and twin similarity (group run) by correlating each PROFUMO mode within each monozygotic twin pair using modes from separate subject runs (#2 above) or subject-specific estimates from the group run (#1 above), respectively. For comparison, we also investigated between-subject similarity (individual runs) and between-subject similarity (group run), reflecting the same across-subject correlations for all possible non-twin pairs of participants. All mode stability measures were performed on each spatial map (correlation across 91,282 vertices), on the temporal connectivity matrix (correlation across the 66 edges in the lower triangle, where edges represent the partial correlation between mode timeseries), and on the spatial overlap matrix (correlation across the 66 edges in the lower triangle, where edges represent the correlation between mode maps; see Section 2.5 and Bijsterbosch et al., 2019).
2.5 Spatial overlap measures
To quantify spatial overlap across all 12 modes, we generated a spatial overlap matrix by estimating the Pearson’s correlation coefficients across grayordinates between all possible pairs across the 12 network maps, as developed in Bijsterbosch et al. (2019). Correlation coefficients were z-transformed prior to averaging across individuals. Stability of individual specific spatial overlap matrices was calculated as described in Section 2.4. Importantly, this spatial overlap matrix definition of network overlap offers a threshold-free estimation of overlapping network organization.
To estimate spatial maps of network overlap in individual participants, each mode was binarized using a threshold of 1, and the number of overlapping modes were counted at each vertex. Each vertex was multiplied by its corresponding cortical area to estimate the spatial overlap area shown in Figure 7B. A threshold of 1 represents the 96.5th percentile across all map weights and therefore offers a relatively conservative estimate of spatial overlap that is driven only by vertices with strong network contributions.
2.6 Focusing on 2-network overlap
To systematically compare our hypotheses regarding network overlap, we focused on spatial overlap between pairs of networks. For each individual and each possible network pair (12*11/2 = 66), we identified vertices uniquely associated with network 1 as those vertices with a weight of 1 or greater for network 1, a weight of less than 0.1 for network 2, and a summed weight across all other 10 modes of less than 0.1 (Fig. 2B). Negative map weights were not included in the estimation of network overlap, because the inclusion of negative vertices would dilute averaged timeseries by canceling out positive vertices. A similar procedure was used to identify vertices uniquely associated with network 2. To locate the spatial overlap region, we identified vertices having a weight of 1 or greater for both network 1 and network 2 and a summed weight across all other 10 modes of less than 0.1. To ensure that findings were not driven by specific thresholding, all nine possible pairings between three thresholds for network vertex weights [0.75 1 1.25] and three thresholds for summed weights [0.02 0.1 0.5] were tested. This procedure was performed using mode maps estimated from the individual subject PROFUMO runs (#2 above). PROFUMO spatial map values in our data ranged from -3.6 to 7.1, and the applied threshold of 1 represents the 96.5th percentile. Hence, thresholds were applied to the PROFUMO spatial maps as estimated by the PROFUMO algorithm without further transforms. As such, we focused specifically on two-mode overlap vertices, by removing vertices with significant contributions of additional modes (Fig. 2A). We focus on 2-network overlap because it offers a tightly controlled test-bed for the hypothesis testing element of our work. Future work may expand into more complex overlapping network organization. Overlap regions were defined based on the data-driven approach described above without any explicit exclusions, and the resulting overlap patterns (see Supplementary Fig. S1) comprehensively cover all areas of potential interest.
The number of network pairs for further investigation was reduced from 66 (all possible pairs) by selecting only those networks pairs in which the spatial overlap regions contained at least 25 vertices for at least half of the participants (n >= 10). This network pair selection was performed to focus on pairs with robust and replicable overlap, and to reduce the computational demands for subsequent analyses. Out of the 66 possible network pairs, 20 pairs were selected for further analysis (see Supplementary Fig. S1).
For each of the 20 network pairs, three summary timeseries were calculated (Fig. 2B) by averaging across: all vertices uniquely associated with network 1 (“N1”), all vertices uniquely associated with network 2 (“N2”), and all vertices in the overlap region (“O”). Each timeseries was standardized to a mean of zero and a standard deviation of 1 within each of the 12 runs. Any participants who did not have any vertices in the overlap region based on these criteria were excluded from subsequent analyses.
2.7 Generating semi-simulated data for the overlap region
Multiple mechanistic hypotheses might in principle explain the occurrence of spatial overlap in stationary maps estimated using PROFUMO, and it is currently unknown what drives the observed (apparent) spatial overlap. To address this issue, we generated semi-simulated versions of overlap timeseries to test different hypotheses as described below. The semi-simulated version of overlap timeseries were compared to the original overlap timeseries using direct timeseries correlation, frequency characteristics (see Section 2.8 for further details), and by fitting a general linear model (GLM; see Section 2.9 for further details).
2.7.1 Network switching hypothesis
One hypothesis is that the overlap region may dynamically switch network alliance between networks 1 and 2 over time within a scanning run (Fig. 1B). Such dynamic switching would appear as overlap when using stationary methods such as PROFUMO, which effectively average across time. To test this hypothesis, we semi-simulated four versions of the overlap timeseries using the real N1 and N2 timeseries for each participant and run:
Switch 50 TRs (Fig. 2D) starts with the first 50 timepoints from N1, then contains the second 50 timepoints from N2, then the third 50 timepoints from N1 and so on. Notably, potential phase shifts between original and semi-simulated data (i.e., mismatches in the order between N1 and N2) may impact the correlation between the original and semi-simulated versions of the timeseries. However, such phase discrepancies would not impact our second comparison measure based on the GLM, because the resulting GLM beta weights will capture order effects.
Switch 25 TRs (Fig. 2D) is the same as above, but switching between N1 and N2 every 25 timepoints.
Switch 10 TRs (Fig. 2D) is the same as above, but switching between N1 and N2 every 10 timepoints.
Max switching (Fig. 2D) assigns each timepoint in the semi-simulated overlap timeseries as the maximum from either N1 or N2 based on whichever datapoint (in real data timeseries N1 and N2) is higher for a given TR.
2.7.2 Coupling hypothesis
Another hypothesis is that the overlap region is integrating data from both network 1 and network 2 at each TR (Fig. 1C). To test the coupling hypothesis, we semi-simulated two versions of the overlap timeseries using the real N1 and N2 timeseries for each participant and run:
2.7.3 Spatial mixture hypothesis
A final hypothesis is that the overlap region is a spatial mixture of vertices linked to network 1 and network 2 (Fig. 1A), akin to the concept of network interdigitation (Braga & Buckner, 2017) or within-region gradients (Haak et al., 2018). To test this hypothesis, we semi-simulated two versions of the overlap timeseries using the real vertex timeseries uniquely associated with either network 1 or network 2. In contrast to the other semi-simulated versions of the overlap timeseries described above, this method does not use the mean N1 and N2 timeseries and instead repeats the averaging across vertices.
Spatial random mixture (Fig. 2C) assigns half of the vertices in the overlap region to a randomly chosen vertex out of those uniquely associated with network 1 (with replacement) and assigns the other half of the vertices in the overlap region to a randomly chosen vertex out of those uniquely associated with network 2 (with replacement). The semi-simulated overlap timeseries is then averaged across all vertices and standardized as above.
Spatial interdigitation (Fig. 2C) assigns each vertex in the overlap region based on whether the spatial weight for that vertex was higher for network 1 or for network 2. If the spatial weight for M-FPN 1 is higher, the timeseries of the vertex uniquely associated with network 1 with the grayordinate that is closest in vectorized indexing is assigned to that overlap vertex (with replacement). The semi-simulated overlap timeseries is then averaged across all overlap vertices and standardized as above. We refer to this option as “spatial interdigitation” because the PROFUMO spatial weights reflect spatially contiguous areas (see Supplementary Fig. S2).
Although the 7 T data were acquired at higher resolution than the 3 T data, we analyzed data on the same 32k mesh (i.e., matched resolution). A deeper investigation into high-resolution 7 T network organization is beyond the scope of this work.
2.8 Frequency characteristics of semi-simulated timeseries
Fourier transforms were performed on the original overlap timeseries and each of the eight semi-simulated overlap timeseries to compare the resulting frequency characteristics. Fourier transforms were performed separately for each participant and each run using only the 3 T runs (8 per individual) for ease of comparison due to matched TR and timeseries length. Resulting power spectra were normalized to a maximum of 1 for ease of comparison.
2.9 Hypothesis testing using the general linear model
The advantage of the semi-simulated timeseries approach above is that it offers a hypothesis-testing framework with direct comparisons across hypotheses based on timeseries correlations between the semi-simulated and original overlap timeseries. The general linear model (GLM) offers a complementary approach based on model fit. Advantages of the GLM approach include data-driven estimation of coefficients. These data-driven parameter estimates enable non-equal contributions from network 1 and network 2 in the linear coupling hypothesis, and can address ordering effects and window-duration challenges in the switching hypotheses. As such, we adopted a GLM framework to corroborate or refute results from the semi-simulated approach. For all GLMs below, the dependent and independent variables were all demeaned and variance normalized prior to model fit. The GLMs were fit separately for each participant and each network pair, and GLMs were also fit separately to data from each of the 12 runs. To enable model comparisons across hypotheses, we adopt adjusted across a range of models described below. Importantly, the adjusted is used to account for differences in model complexity.
2.9.1 Network switching hypothesis
To test the switching hypothesis, two separate GLMs were fit to each windowed segment, using the following equations:
Windows with Network 1 alliance:
Windows with Network 2 alliance:
Here, refers to the overlap timeseries, refers to the timeseries from network 1, and refers to the timeseries from network 2 (all derived from the original data, averaged across vertices). indicates the window, which is modeled as a sliding window of 10 timepoints shifted in steps of 1 timepoint. Note that there is no need to model longer sliding window durations, because any switching behavior in longer windows will be captured across multiple shorter windows. For each window, the adjusted of the best performing model (either network 1 alliance or network 2 alliance) was recorded, and the resulting vector of was averaged across windows. If the overlap region is dynamically switching alliance over time, then any individual segment should be well explained by either N1 or N2, and the resulting (selected as the highest between the two models above) should therefore be equal to or higher to for other GLMs described below.
2.9.2 Coupling hypothesis
To test the coupling hypotheses, separate GLMs were fit to compare the linear, combined (including linear and nonlinear terms), and nonlinear coupling hypotheses. Models were fit to the overall timeseries using the equations:
Linear hypothesis:
Combined linear and nonlinear model: +
Nonlinear hypothesis (interaction only):
The nonlinear GLM models the interaction term, which is calculated by point-wise multiplying N1 by N2, after setting the minimum of both N1 and N2 to zero. The for each model was recorded and compared across other hypotheses.
2.9.3 Spatial mixture hypothesis
To test the spatial mixture hypotheses, separate GLMs were fit to each vertex timeseries in the region of network overlap using the equation:
Spatial mixture hypothesis:
Here, the indexing of on refers to the vertex under investigation. To enable direct comparison with other hypotheses, the resulting was averaged across all . The resulting spatial patterns were not further investigated due to poor performance (see Results, Section 3.6).
3 Results
3.1 PROFUMO mode maps
Out of the total of 20 modes, 12 met the test-retest and subject-group criteria to be considered for further analyses (see Section 2.3). The 12 modes (Fig. 3) covered well-known occipital (visual), pericentral (somatomotor), dorsal frontoparietal (attention), lateral frontoparietal (control), and medial frontoparietal (default) networks (Uddin et al., 2019). These 12 modes were highly replicable across different dimensionalities of PROFUMO (Supplementary Fig. S3). Individual participants had on average 1.4 missing modes (range 0-4; Supplementary Fig. S4). Group maps were stable across the group-PROFUMO using all data and group-PROFUMO using only 12 runs to match subject analyses, albeit with lower weights in the 12-run results reflecting the reduction in SNR (Fig. 4 top row). Mode maps for two example participants derived from single-subject PROFUMO reveal detailed individual specific organization that closely matches maps from the same participants derived from the classic PROFUMO group analysis (Fig. 4 middle and bottom rows). The example participants were chosen as non-twin individuals with a complete set of 12 modes and are representative in terms of all other indices (as shown by highlighting the blue and red data points in Figs. 4, 5, 7 & Supplementary Figs. S4, S6, S7, S8). Subject-specific networks and twin comparisons are shown in Supplementary Figure S5.
3.2 PROFUMO mode stability
An initial objective was to determine whether PROFUMO can reliably estimate resting state maps using only data from a single participant. Figure 5 shows that the highest similarity occurs for single-subject PROFUMO spatial maps (test-retest reliability and within-subject similarity in columns 1&2, respectively; mean r = 0.80 ± 0.18); the next highest similarity occurs for twin correlations (Fig. 5 columns 3&4; mean r = 0.70 ± 0.14); and the lowest similarity occurs between non-twin participants (Fig. 5 columns 5&6; mean r = 0.59 ± 0.15). Each similarity pattern shows a “tail” of networks having lower stability (Fig. 5), which on further inspection appeared to be distributed across participants and across networks. When all modes are included (i.e., not only the 12 selected modes and not removing “missing” individual modes), these tails are further expanded but the findings described above regarding comparisons between columns remain evident (see Supplementary Fig. S6). Single-subject PROFUMO spatial maps derived using the 12 runs for the participant (i.e., not including data from other participants) were also highly similar to the estimated subject maps derived from the group data informed by all 20 subjects with 12 runs each (Fig. 5 column 2; mean group-subject r = 0.86 ± 0.11), suggesting reasonable correspondence at least at this dimensionality. Spatial maps estimates informed by group data achieved slightly higher similarity compared to maps estimated from individual runs (Fig. 5 column 4 > 3 and column 6 > 5), revealing the impact of the group prior when performing classic hierarchical PROFUMO. Stability estimates for second-order statistics, including temporal connectivity matrices and spatial overlap matrices, are shown in Supplementary Figures S7 and S8.
3.3 Amount of data needed to estimate single-subject modes using PROFUMO
To test how much data from an individual participant were needed to obtain good estimates from single-subject PROFUMO, we systematically varied the number of timepoints (in increments of 1/12th of the run timepoints), and compared the resulting spatial maps to the results that included all timepoints. Mean similarity across the selected (non-missing) mode maps increased when the number of timepoints were increased (Fig. 6) and was near asymptotic above approximately 5000 TRs (approximately 1 hour of data per participant).
3.4 Spatial overlap
Consistent with previous work (Bijsterbosch et al., 2019), our results indicate substantial spatial correlation across PROFUMO modes (Fig. 7A). Figure 7A shows the spatial overlap matrix calculated by correlating pairs of weighted PROFUMO group spatial maps. Notably, spatial overlap is ignored in the find-the-biggest overview in Supplementary Figure S5, in which vertices are assigned to a single PROFUMO mode with the highest spatial weight to enable concise visualization of all subject maps. Spatial overlap was primarily localized in the lateral parietal cortex and posterior cingulate - precuneus regions (Fig. 7C, D), consistent with prior work (Bijsterbosch et al., 2019). Spatial overlap maps for two example participants are shown in Figure 7C and D; Supplementary Figure S9 shows maps for all participants and across twin pairs, and Supplementary Figure S1 shows maps for the selected set of 20 2-network overlap pairs. The 2-network spatial area estimates were not associated with DVARS (the temporal derivative of the root mean square variance (Power et al., 2012); r = -0.06, puncorrected = 0.81), suggesting that spatial overlap was not sensitive to head motion. Conversely, a trend-level association was observed between 5-network spatial area and DVARS (r = -0.45, puncorrected = 0.05), although we note that our sample size is insufficient to reliably model individual difference (Marek et al., 2022).
3.5 Semi-simulated data for 2-network overlap
We derived semi-simulated versions of the overlap timeseries based on combinations of network 1 and network 2 timeseries (Fig. 2C–E). We assessed how similar the semi-simulated timeseries for the overlap region were, compared to the original overlap timeseries (O) by estimating the Pearson’s correlation coefficient separately for each network pair, each subject, and each run. The highest similarity was observed for the linear additive coupling hypothesis, which achieved a median correlation of 0.783 (Fig. 8). This result was significantly different from the next highest correlation (median of 0.767), observed for the nonlinear multiplicative coupling hypothesis (T = 6.3, p = 3.2*10-10, estimated after z-transformation of the correlation values). These findings were observed reliably across all thresholds used for network inclusion and exclusion (see Supplementary Fig. S10).
Across the 20 different network pairs, the linear additive semi-simulated timeseries achieved the highest correlation for 16 network pairs, and the nonlinear multiplicative semi-simulated timeseries achieved the highest correlation for the remaining 4 network pairs (see Supplementary Fig. S11). Notably, each of the 4 network pairs with highest correlations for the nonlinear multiplicative semi-simulated timeseries involved the L-FPN 1 (paired with M-FPN 2, L-FPN 2, L-FPN 3, and L-FPN 4 respectively). The L-FPN 1 network includes a specific subregion (“POS2”) of the parietal-occipital sulcus that is highly distinctive from neighboring regions in its myelin, thickness, connectivity, and task activation (Glasser et al., 2016). As such, the L-FPN 1 network (and POS2 area in particular) warrants future research into their distinctive features, including nonlinear origins of network overlap.
Frequency characteristics were highly similar between the original overlap timeseries and all semi-simulated timeseries (see Supplementary Fig. S12), although spatial interdigitation timeseries exhibited a slightly raised tail and the temporal switching hypothesis (every 10 TRs) has somewhat increased power at 0.1 Hz.
3.6 General linear model results
The results from the semi-simulated timeseries analysis revealed that the linear coupling hypothesis was most strongly supported, followed by the nonlinear coupling hypothesis. We subsequently adopted a GLM-based approach to corroborate these findings. Results confirm that the linear coupling hypothesis achieved the best model fit (median = 0.65), which was significantly higher than the nonlinear coupling hypothesis (interaction only; median = 0.59, T = 13.2, p = 4.5*10-39; Fig. 9), and was not significantly lower than the combined model with the nonlinear term added (median = 0.65; T = -0.23, p = 0.82).
4 Discussion
Our first aim was to determine whether weighted resting-state spatial networks can be robustly derived from single-subject data without sacrificing correspondence. The results revealed 12 resting-state networks with high test-retest reliability and remarkably good incidental (i.e., non-modeled) correspondence with group-informed network estimates (Fig. 4). Notably, the test-retest reliability of these individual-specific network maps (cosine similarity approximately 0.8 based on Fig. 10A in Harrison et al., 2020) matches or exceeds many previous estimates of reliability of resting-state metrics (e.g., average intraclass correlation 0.29 for temporal functional connectivity from meta analysis across 25 studies; Noble et al., 2019) (Andellini et al., 2015; Braun et al., 2012; Dutt et al., 2022; Fiecas et al., 2013; J. Wang et al., 2017; J.-H. Wang et al., 2011; Liao et al., 2013; Nemani & Lowe, 2021; Noble et al., 2019; Termenon et al., 2016; Yang et al., 2021). These findings support the possibility of future personalized psychiatry approaches where data from an individual would be separately analyzed and compared to a reference cohort to inform clinical decision making. It is possible that increasing the dimensionality to extract more finer-grained resting-state networks (beyond the 12 selected networks out of 20 used here) may reduce test-retest reliability and correspondence. However, we showed high stability of the selected networks across higher dimensionalities of PROFUMO (Supplementary Fig. S1), and previous work in networks defined with independent component analysis reported good reliability up to at least 150 networks (Ma & MacDonald, 2021). Importantly, the functional organization of the brain can meaningfully be studied at multiple levels of complexity along its organizational hierarchy (Bijsterbosch et al., 2021). Here, we specifically chose to investigate the organization of macroscale networks (i.e., low dimensionality) because it describes functional organization in relation to widely studied canonical networks that are consistently observed across datasets, analysis methods, and states (Smith et al., 2009; Uddin et al., 2019; Yeo et al., 2011).
Our second aim was to determine the degree of spatially overlapping network organization. Our findings confirm the presence of extensive and stable network overlap in networks estimated from single-subject data. Patterns of spatially overlapping network organization vary extensively across individuals (Fig. 7 & Supplementary Fig. S6), and previous work has shown that individual differences in these overlap patterns are strongly associated with behavior (Bijsterbosch et al., 2019). As such, the lower dimensional spatial overlap matrix (Fig. 7A) may provide a key summary measure of behaviorally-relevant aspects of spatial organization, while reducing the multiple-comparison burden of vertex-wise analysis of spatial organization.
Spatial overlap between individual-specific weighted resting-state network maps also offers a complementary approach to investigate hub regions. As opposed to traditional hub identification methods that rely on temporal correlations, our weighted network approach emphasizes the role of shared brain regions as part of the spatial organization of functional networks. It is currently unclear whether these distinct temporal versus spatial definitions of hub regions have dissociable or shared neurobiological implications. Some evidence suggests that spatial networks and hubs may maintain functional integrity over time, potentially for homeostatic purposes (Laumann & Snyder, 2021), whereas flexible hubs that can rapidly change their temporal connectivity may provide coordination and switching functions to support cognitive control (Cocuzza et al., 2020; Cole et al., 2013; Gordon et al., 2018). Consistent with this hypothesis, we have shown that spatial overlap was more strongly associated with stable trait-like behavior (Bijsterbosch et al., 2018), whereas temporal correlations tracked transitions between sensorimotor states (Harrison et al., 2020). However, our prior work also suggests that effects of spatial network organization and overlap are observed as temporal correlation estimates when unaccounted for in the parcellation (Bijsterbosch et al., 2018, 2019), which may indicate that the distinction between spatial and temporal hubs could be purely analytical, driven by differences in model cost functions and priors, without distinct neurobiological interpretations. Hence, spatial network overlap and temporally strongly connected nodes may serve distinct functional purposes (e.g., homeostasis versus switching), or may represent alternative analytical estimates of the same underlying neural phenomenon. Further work is needed to gain insight into spatial versus temporal hub definitions and their neurobiological functions.
Our third aim was to systematically test different hypotheses (spatial mixture, dynamic switching, and coupling hypotheses; Fig. 1) regarding the nature of brain regions in which multiple resting-state networks appear to overlap. Our findings supported the linear coupling hypothesis for the 2-network overlap case (Figs. 8 & 9). Future work investigating more complex patterns of overlap across three or more networks might provide support for a different hypothesis. Given that regions of network overlap appear to be actively engaged in processing information from all contributing networks, this suggests that network overlap may play an important role in the integration of information across multiple brain systems. Notably, the results did not conclusively differentiate between linear and nonlinear versions of the coupling hypothesis. This differentiation will be of interest for future research because the linear additive hypothesis may suggest that multiple networks may coexist within overlap regions without influencing one another. Specifically, the combined signal may indicate that both networks behave as they do in non-overlap regions without any cross-network modulation or integration. The nonlinear hypothesis, on the other hand, requires cross-network integration. Prior work has reported that nonlinear binding between multiple task conditions in conjunction hubs (i.e., brain regions that selectively integrate activations) was essential for predicting task activation patterns from functional connectivity data (Ito et al., 2022). Intriguingly, performance for the nonlinear multiplicative hypothesis was not far behind the additive linear hypothesis, which may indicate the presence of network interaction. Although the biological mechanisms of network overlap and additive/multiplicative coupling are unclear, we hypothesize that regions of network overlap do not constitute stand-alone cortical areas but instead play an important role in cross-network integration, for example through the synchronization of local field potentials in band-limited frequency ranges or their enveloppes (Acland et al., 2023).
Despite offering novel insights into weighted resting-state networks estimated within individual participants, this work also has several limitations. First, we investigated network organization at a PROFUMO dimensionality of 20 and investigated 12 resulting networks of high stability. Notably, the 12 networks under investigations were obtained with very high stability across a range of PROFUMO dimensionalities (Supplementary Fig. S1). Nevertheless, network overlap may behave differently at different dimensions of network decomposition, and both higher (Farahibozorg et al., 2021) and lower (Lenzini et al., 2023) network dimensions are of interest for future research. Second, our vertex assignments to network 1, network 2, and overlap regions involved thresholding, which is necessarily a simplification of the weighted network organization. However, the strict threshold for exclusion and inclusion of vertices enforced a relatively conservative definition of overlap that was necessary to make the results interpretable. Importantly, we recommend the use of the unthresholded spatial overlap matrix (Fig. 7A) for brain-behavior investigations of network overlap. Third, the estimation of network overlap, whether based on weighted or thresholded maps, interacts with the contrast-to-noise ratio (CNR) of the data, which varies across the cortex and is higher in the lateral parietal region. As such, potential future efforts to develop a conclusive map of network overlap in the brain should perform careful scaling/normalization to remove any bias resulting from CNR variability across the cortex. Furthermore, individual difference studies should match the amount of data used across individuals to avoid biased estimates of network overlap that may arise from CNR variability across individuals. Fourth, our findings indicating support for the coupling hypothesis are consistent with the underlying outer product model used in PROFUMO (and similar methods such as ICA). In future work, it will be of interest to test mechanistic hypotheses of network overlap using alternative approaches to define the overlap region that do not rely on the outer product model. Notably, alternative tools for the definition of network overlap (such as Karahanoğlu & Van De Ville, 2015; Li et al., 2017) might render different overlap results and might have different reliability and data-needs. Fifth, it remains a possibility that spatial mixing and/or dynamic switching may occur at finer spatial and temporal resolutions that cannot be resolved using resting-state MRI data. Furthermore, it is possible that temporally lagged correlation structure may play a role in regions of spatial overlap, which can be challenging to accurately discern using functional MRI data (Smith et al., 2011) but may be feasible in some situations using deconvolution approaches (Mill et al., 2017). Although beyond the scope of the current paper, we plan to test spatial mixing, dynamic switching, and temporally lagged hypotheses at sub-MRI scales in future research using invasive recording techniques in non-human primates.
Taken together, we showed that weighted resting-state networks derived from single-subject data are stable, correspond closely to group-informed networks, and capture overlapping network organization, and are therefore important targets for clinical biomarker research. We also showed that overlapping network organization is indicative of coupling between networks, providing a mechanistic hypothesis for the functional role of these regions.
Ethics
Secondary data analysis of HCP data for this project was covered under IRB 202105086.
Data and Code Availability
The data used in this paper are freely available from https://db.humanconnectome.org/ (requires free registration). All code for this paper is available on GitHub: https://github.com/JanineBijsterbosch/Individual_PROFUMO. Brain maps from Figures 3, 4, and 7 can be viewed in Balsa (including maps for all participants): https://balsa.wustl.edu/study/gm40X.
Author Contributions
J.D.B. conceived of the study, performed all analyses, and wrote the manuscript. S.-R.F. provided support regarding the PROFUMO algorithm. M.F.G. and D.V.E. helped with HCP aspects and interpretation. L.H.S. contributed to the hypothesis generation and funding. M.W.W. and S.M.S. helped with the analysis design and interpretation. All authors reviewed the manuscript.
Funding
This project was supported by the NIH (NINDS R34 NS118618; co-PIs J.D.B. & L.H.S.). M.F.G. and D.V.E. are supported by NIMH R01 MH60974. S.M.S. and M.W.W. are supported by Wellcome Trust (215573/Z/19/Z). S.R.F. is supported by the Royal Academy of Engineering under the Research Fellowship programme (RF2122-21-310).
Declaration of Competing Interest
The authors declare no competing interests.
Supplementray Materials
Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00046