Revealing dynamical organization of the brain

To test the efficacy of our approach in estimating a representation of the brain’s dynamical organization and in capturing transitions in the whole-brain activity, we employed already collected fMRI data from Gonzalez-Castillo et al.34. These data were collected while participants performed in a Continuous Multitask Paradigm (CMP) with experimentally constrained transitions from one task to the next.

After standard fMRI preprocessing, each participant’s 4D fMRI data were first transformed into a matrix, such that the rows corresponded to the individual time frames (or volumes) and the columns corresponded to the intensity value at each voxel. Thus, each row of this matrix represents the entire brain volume at any time point during the session (Fig. 1a). The TDA-based Mapper is next employed on this matrix to perform four steps—filtering, binning, partial clustering, and finally constructing the shape graph (Fig. 1b−e). Although data-driven optimization was employed to find the best set of parameters for each of the Mapper steps, we observed that the presented results are stable in the face of extensive parameter perturbations (see Reliability of shape graphs).

Fig. 1 Application of Mapper to 4D fMRI data. a Pre-processed four-dimensional fMRI data from each participant was fed into the analysis. For each participant, the entire data matrix (i.e., #TRs × #Voxels) was used for analysis. b Non-linear dimensionality reduction was done during the filtering step to project fMRI data into a lower two-dimension set (represented by dimensions (f, g)). c Two-dimensional binning was then performed by dividing the lower-dimension space into smaller bins (determined by the resolution parameter, R) with certain overlap (determined by the gain parameter, G). d Partial clustering was then performed to get a compressed representation by collapsing data into fewer nodes, where each node represents a cluster, and the size of each node depicts the number of data points inside each cluster. e After clustering, nodes that share data points (i.e., time frames in this case) are linked together with an edge to create the final compressed combinatorial representation (or graph) Full size image

As a first Mapper step, we applied a filtering (or dimensionality reduction). This step is similar to the standard dimensionality reduction techniques used in the machine learning literature. However, unlike traditional linear dimensionality reduction techniques, like principal component analysis (PCA) or multi-dimensional scaling (MDS), we employed a nonlinear dimensionality reduction method using a variant of stochastic neighborhood estimation (SNE36,37). Nonlinear methods like SNE allows for preservation of the local structure in the original high-dimensional space after projection into the low-dimensional space, which is typically not possible with linear methods like PCA or MDS36. Thus, the time frames with similar activation patterns in the original high-dimensional space will be projected closer to each other in the reduced dimensional space (Fig. 1b).

To encapsulate the low-dimensional representation generated by the filtering step, Mapper employs binning (or partitioning) (Fig. 1c), followed by partial clustering within each bin. The binning step partitions the low-dimensional space into overlapping bins by using two parameters—number of bins (or resolution (R)) and percentage of overlap between bins (or gain (G)). Within each bin, single-linkage clustering is performed to condense the time frames into a set of one or more clusters (Fig. 1d). This step results in a compressed representation, as fewer points (or clusters) are now required to represent the data as compared to the entire set of time frames. On average, for each participant, the compressed representation contained ~279 points (SD = 60) (as opposed to 1017 acquired time frames). It is important to note that this clustering step is different from traditional temporal smoothing as the time frames within each cluster are not averaged and the mapping between individual clusters and their corresponding time frames is preserved.

Finally, to generate a combinatorial object or shape graph from the low-dimensional compressed representation, the Mapper treats each cluster as a node in the graph and connects these nodes with an edge if they share time frames (Fig. 1e). The final shape graph can be conceptualized as a low-dimensional depiction of how the brain dynamically evolved across different functional configurations during the scan. While the actual interpretation of the latent variables associated with the projected low-dimensional space may differ across participants, the topological relationships encoded by the shape graph itself are interpretable and comparable across participants (as we shall discuss later).

To reveal the underlying temporal structure in the CMP dataset, the aforementioned Mapper approach was applied to each participant in the CMP dataset. For qualitative analysis, we annotated the nodes in these shape graphs with colors based on the corresponding task at each time frame (Fig. 2a). Further, if a node contained time frames from multiple tasks, we visualized that node using a pie chart denoting the proportion of time frames that belong to each task within such node (Fig. 2b). Graph theoretical metrics were next used to quantify the topological properties of each participant’s shape graph.

Fig. 2 Revealing the shape of brain’s dynamical organization. a Depicts the shape graph for one of the representative participants (S 01 ) from the CMP dataset. The shape graph was annotated using colors and pie-chart visualization scheme to depict how the tasks were represented in each shape graph. b Shows a zoomed-in version of the densely connected region of the shape graph to show the use of pie-chart visualization Full size image

Quantifying the mesoscale structure of shape graphs

Graph theory (or Network Science) is currently widely used in the field of neuroscience to provide summary statistics of the complex interactions between different entities or nodes. While interesting insights can be captured by analyzing properties of each node or edge in the network (i.e., at the local scale) or by analyzing the network as whole (i.e., at the global scale), the intermediate (or mesoscale) properties appear particularly well suited for analyzing and comparing the structure of complex networks38,39. In particular, considerable effort has gone into identifying two distinct types of mesoscale structures in a variety of complex networks. The first and perhaps the most widely used mesoscale structure is the community structure, where cohesive groups called communities consist of nodes that are densely connected to other nodes within communities while being only sparsely connected to nodes between communities40. In the context of shape graphs, the presence of communities could represent a modular organization with specialized whole-brain functional configurations for different types of information processing (or tasks). An increasingly second most typical mesoscale structure is the core−periphery structure41. Here, one attempts to determine the core nodes, which are not only densely connected to each other but are also central to the entire network. A presence of core nodes in shape graphs could indicate whole-brain functional configurations that consistently occur throughout the scan. For example, core nodes could represent neural processes related to task-switching that the brain consistently passes through during a multitask experimental paradigm. The peripheral nodes on the other hand are only sparsely connected. The examination of the core–periphery structure of a graph could reveal the overall arrangement of the network39. It is important to note that in the real world, networks can have both communities and core−periphery structures and hence it is desirable to investigate both simultaneously.

To exemplify how topological properties of shape graphs can provide behaviorally relevant information at the single-participant level, we estimated both community and core−periphery mesoscale structures. Briefly, to estimate the degree or quality of the community structure in shape graphs, we assessed the widely used quality function Q mod 42. The community assignment for each node in a shape graph was chosen to be one of the four tasks (i.e., Rest, Memory, Video, and Math) based on the majority of time frames contained in the node that belonged to the respective task. Across the CMP dataset, we observed participants’ shape graphs with varying degree of modularity (ranging from Q mod = 0.37 to 0.61 with a mean = 0.48 and SD = 0.07; Fig. 3a). Remarkably, the degree of modularity was observed to be associated with task performance across the three CMP tasks (%correct r = 0.56, p = 0.016; Fig. 3b), such that high modularity was associated with better performance on the CMP tasks. Thus, highlighting that participants with a higher degree of community structure in their shape graph better performed across different cognitive tasks during the CMP. In other words, participants with specialized whole-brain configurations for different tasks were those with the highest overall task performance. We examined this claim using an independent validation analysis (see Anchoring topology of shape graphs into anatomy). Please note that because the quantification of community structure was done on the shape graph as a whole (i.e., across tasks), we combined task performance measures across the working memory, video and math tasks.

Fig. 3 Quantifying the mesoscale community structure of shape graphs. a Shows two shape graphs from two representative participants (S 14 and S 07 ). As depicted, S 14 had a low modularity score, i.e., nodes with different task types are connected to each other without any preference for same task nodes. On the other hand, shape graph from S 07 depicts high modularity structure with nodes preferentially connected to other nodes of same task type. b The modularity score (Q mod ) was observed to be associated with task performance, such that higher modularity in the shape graph was associated with better task performance Full size image

To quantify the core−periphery structure in each participant’s shape graph, we employed the generalized Borgatti and Everett41 algorithm that provides a coreness score (CS) for each node. This algorithm assigns CS along a continuous spectrum with nodes that lie most deeply in a network core with a CS ~1 to those that are in the periphery with a CS ~039. Figure 4a presents shape graphs annotated (or colored) by the task type as well as CS for two representative participants. Remarkably, across all participants, the nodes containing resting state time frames were mostly represented in the peripheries (mean CSRest [SD] = 0.15 [0.06]), while the nodes containing time frames from cognitively demanding tasks mainly lied relatively deeper inside the shape graph (mean CSW.M. [SD] = 0.28 [0.04]; mean CSMath [SD] = 0.30 [0.04]; and mean CSVideo [SD] = 0.22 [0.06]). One-way ANOVA revealed a significant effect of the task (F(3,51) = 24.06, p < 0.0001), such that CSRest was observed to be significantly lower than the CS of other three tasks, while coreness scores for the working memory and math tasks were similar but higher than that of the video task. This result indicates more consistency in the whole-brain functional configurations was present during math and memory task as compared to the less demanding resting state.

Fig. 4 Quantifying the mesoscale core−periphery structure of shape graphs. a Shows shape graphs for two representative participants (S 14 and S 01 ). In the left column, we annotated the nodes using task-type and in the right column the nodes are annotated using the coreness scores (CS). As evident, the more central and densely connected nodes have high CS and more sparsely and peripheral nodes have less CS. b Shows boxplots of CS derived from the CMP data as well as from the three null models (neuRosim, phase randomization with random, and constant phase). Across the four CMP tasks, nodes containing resting time frames had the lowest CS, while the nodes containing time frames from cognitively demanding tasks mainly lied relatively deeper inside the shape graph (i.e. higher CS). The CS of nodes with resting time frames was also lower than all the three corresponding null models (ps < 0.001), and the CS of nodes with working memory or math frames was higher than all the three corresponding null models (ps < 0.005). No significant difference was observed for the CS of nodes with time frames from the video task and the corresponding null models. c Shows the brain regions that were positively associated with coreness scores, i.e., higher the CS, higher the activation in these regions. The group-maps (in red-yellow scale) for positive relation with CS are overlaid on the task-associated meta-analysis maps from NeuroSynth44. As shown, higher CS was associated with higher engagement of task-related brain regions Full size image

To test the validity of the observed non-trivial arrangement of resting state nodes in the peripheries while the cognitively demanding nodes in the core, we employed three different null models. Overall, the CS of nodes with resting time frames was observed to be lower than all the three corresponding null models (ps < 0.001), while the CS of nodes with working memory or math frames was higher than all the three corresponding null models (ps < 0.005). No significant difference was observed for the CS of nodes with time frames from the video task and the corresponding null models (Fig. 4b).

Anchoring topology of shape graphs in anatomy

To ground the shape graphs and their properties into neurophysiology, we provide three approaches that attempt to reveal the underlying patterns of brain activity putatively responsible for the observed topological features. In the first approach, we use spatial mixture modeling (SMM)43 to reveal changes in brain activation maps from one time frame to the next. The SMM approach includes fitting a mixture of distributions and using a spatial Markov random field to regularize the labeling of voxels into null, activated or deactivated. Thus, for each node in the shape graph and the containing time frames, we generated whole-brain activation (and deactivation) maps. To interactively examine the temporal variations in these activation maps, we developed a web-tool (Supplementary Figs. 2−3 and Supplementary Movie 1). The web-tool allows for a better interpretation of the shape graph and its topological properties (i.e., edges and nodes). For example, in real time, a user can move the Time-Frame slider (across time frames) to simultaneously highlight respective nodes in the shape graph; see transitions in corresponding whole-brain activation maps; and observe correlations of the activation maps with known large-scale brain networks9 (Supplementary Movie 1). Thus, allowing an inspection of neurophysiology at the whole-brain level and the highest temporal resolution (limited only by acquisition rate). When using this tool on the CMP dataset, we can observe how nodes associated with the memory task correspond to whole-brain activation maps with significant activity in the dorsolateral prefrontal cortex, and visual cortex; while for the math task we see the engagement of parietal regions previously associated with arithmetic processing (Supplementary Fig. 4). These results validate the cognitive relevance of the activity maps generated from the shape graphs as they identify regions commonly activated by these tasks.

In the second approach, to anchor the overall topological properties of the shape graph into neurophysiology, we utilized the traditional group-based generalized linear model (GLM) analysis. Specifically, we examined the neurophysiological basis for the observed non-trivial mesoscale structure of core–periphery in the shape graphs. For this analysis, the coreness score of each node was mapped back to the individual time frames contained in that node. Thus, if a node has a CS of 0.5, then the time frames contained in that node also received a CS of 0.5. Using multiple regression, the CS for each time frame was entered for each task creating four explanatory variables. Two contrasts were run to examine brain regions that show positive as well as negative association with the coreness scores. The cluster-corrected (Z > 2.3 and FWER p < 0.05) group-level results are shown in Fig. 4c and Supplementary Table 1. For the positive association contrast, during the working memory task, higher coreness scores were associated with increased engagement of the bilateral dorsolateral prefrontal cortex (DLPFC), bilateral insula and lateral occipital cortex, and paracingulate gyrus. Higher coreness scores during the math task were associated with increased engagement of the R. angular gyrus, inferior parietal sulcus areas and the paracingulate gyrus. Lastly, for the video task, higher coreness scores were positively associated with activation in the bilateral fusiform gyrus and right frontal pole. Qualitatively, the brain regions associated positively with coreness scores largely overlapped with regions previously shown to be recruited for the respective tasks (Fig. 4c). This qualitative assessment was performed by overlaying the observed results on the meta-analysis maps generated by NeuroSynth.org44 for each task term (i.e., “working memory,” “arithmetic,” and “object recognition,” respectively). For the negative association contrast, across all three tasks, significant clusters were observed in the posterior cingulate cortex (PCC) and medial prefrontal cortex (Supplementary Fig. 5), suggesting nodes with lower coreness scores (i.e., periphery nodes) were associated with increased activation in the PCC irrespective of the task type. No significant cluster was observed for positive or negative association with coreness scores in the resting state task. Overall, these results suggest that core nodes in the shape graph represent task-related activation and putatively associated cognitive effort, whereas sparsely and peripherally connected nodes in the shape graph represent task-unrelated activation presumably related to task-negative default mode regions.

Lastly, an independent whole-brain FC analysis was run to validate the TDA-derived finding that individuals with higher task performance putatively evoked task-specific brain configurations as compared to individuals with lower task performance. Two groups were formed using a median split based on the overall task performance (low (n = 9) and high performers (n = 9) with a median split at 86.9% accuracy). The whole-brain FC was estimated by sampling data from a set of 264 brain regions to make inferences at the regional and systems level (Supplementary Fig. 12). We used an independent and well-established brain parcellation scheme that was previously identified using a combination of resting-state FC parcellation as well as task neuroimaging meta-analysis45. The FC between each pair of brain regions was estimated using Pearson correlations (r), and was converted using Fisher’s z-transform for further analysis. After estimating whole-brain FC matrices for each task-block, we calculated the similarity between FC matrices derived from different task blocks within each participant. To perform a group-level comparison, the average similarity between FC matrices across tasks, for each participant, was estimated and a two-sample ttest was run to compare the low-performance group from the high-performance group. A significant t-statistic (t(16) = 2.50, p = 0.0236) was observed, such that participants in the low-performance group (compared to the high-performance group) had higher average similarity between FC matrices across tasks (Supplementary Fig. 12B–C). Thus, validating the TDA-derived prediction that participants who performed better across tasks evoked task-specific brain configurations.

Capturing temporal transitions in brain dynamics

To estimate the temporal transitions in the whole-brain activity maps, we converted the TDA-generated shape graph into an adjacency matrix in the temporal domain (i.e., a temporal connectivity matrix (TCM)). It is important to note that a TCM is representing similarity (or connectivity) in time and not in space like the standard FC matrices that represent anatomical region-by-region connectivity. Here, the time frames are considered connected if they share a node or if the nodes containing these time frames are connected by an edge in the shape graph. As evident from Fig. 5a, for a representative participant, the TCM was observed to be modularly organized, with densely connected frames within each task block and across blocks of the same task.

Fig. 5 Capturing temporal transitions at the level individual time frames. a Shows a temporal connectivity matrix (TCM) for participant S 01 from the CMP dataset. A TCM (size: #time frames × #time frames) shows how each time frame is connected (or similar) to all other time frames. The time frames are connected (i.e. non-zero value in the TCM) if they share a node in the shape graph or if the nodes containing these time frames are connected by an edge in the shape graph. We have manually highlighted the task blocks on the TCM to show the task-related modular structure. b Shows that the degree of TCM nodes can capture the transition between tasks. Here we show average degree across all 18 participants as a solid-line, shaded region shows the standard error around the mean. The dotted lines denote detected change points, which mark the transitions for the eight task blocks. c Shows the effectiveness of capturing temporal transitions extracted from the CMP data as compared to transitions extracted from the phase-randomization null model. The blue dashed line denotes onset of task block and the blue dotted line indicates offset of task block. As compared to the null model data (denoted using diamonds), a time frame by time frame analysis revealed that the degree of TCM generated using real data (denoted using + signs) could capture both the onset and offset of each task within a matter of a few time frames. The red inverted triangles on top denote point-by-point significant difference between transitions from real as compared to null data Full size image

Remarkably, by directly estimating the degree (or the total number of connections) at each time frame in the TCM, we could capture the transitions between and within task types at the level of a few time frames. Inherently, a higher degree at any time frame implies greater similarity of that frame with other frames. Thus, during the task blocks, other than resting state block, the evoked activity associated with the stimuli/task would cause the time frames to be highly coherent or similar within each block and across the repetition of the same task (and hence more connected in the TCM), thereby leading to a higher degree value. During the resting state blocks (as well as during between-task instruction periods) the brain activation patterns were driven by intrinsic (and not evoked) activity, which would lead to less coherent or dissimilar patterns and hence a lower degree value. Thus, a task-switch from an evoked task to an instruction period or vice versa would lead to a change in degree values at the level of a few time frames (Fig. 5b, c). As shown in Fig. 5b, the average normalized degree of TCM reveals the between task transitions. Using a standard change point detection algorithm46, we were able to retrieve eight transitions in the mean normalized degree, corresponding to the eight task blocks. Additionally, for each task type, we quantified how quickly the degree of TCM can capture both onsets and offsets of each task block. As compared to the null model generated using a phase randomization technique (see Validation of the shape graphs against null models), a time frame by time frame analysis revealed that the degree of TCM generated using real data could capture both the onset and offset of tasks within a matter of a few time frames (Fig. 5c).

Similar to the coreness score analysis, a one-way ANOVA revealed a significant effect of the task for the degree at each time frame in the TCM (F(4,68) = 104.27, p < 0.0001). The working memory task had the highest degree (mean \({{d}}_{{\mathrm{TCM}}}^{{\mathrm{WM}}}\)=0.71 [0.11]), followed by math (mean \({{d}}_{{\mathrm{TCM}}}^{{\mathrm{Math}}}\)=0.60 [0.13]), then video (mean \({{d}}_{{\mathrm{TCM}}}^{{\mathrm{Video}}}\)=0.31 [0.17]) and lastly, resting state (mean \({{d}}_{{\mathrm{TCM}}}^{{\mathrm{Rest}}}\)=0.15 [0.06]).

Reliability of shape graphs

We performed three reliability analyses. First, for each participant, we split the data into two halves and ran our approach on each half independently. We hypothesized that if our approach is reliable, we should observe similar shape graphs and their properties for the two halves. As hypothesized, the independently generated shape graphs for both halves had similar mesoscale structural as well as temporal properties as observed with the full dataset. In the shape graph from each half of the data, we found: (1) high modularity associated with better performance across tasks; (2) core−periphery structure, such that coreness scores for resting state was lower than cognitively demanding tasks; and (3) between-task temporal transitions at the level of few time frames (Supplementary Fig. 6). These findings suggest that our approach is replicable even when only a part of the data is used to estimate shape graphs.

As a second attempt to estimate reliability, we replicated our approach on an entirely different dataset with distinctive data acquisition parameters, preprocessing steps and experimental design. We used fMRI dataset from 38 unrelated participants from the Human Connectome Project (HCP)34 who performed a working memory task. The working memory task consisted of two runs (~5 min each) and three conditions (2-back, 0-back, and fixation). As compared to the CMP dataset, the HCP dataset had only one task for each session. Thus, instead of examining dynamical organization across tasks, we examined the organization within a single task (but across conditions). Using similar Mapper parameters as used in the CMP dataset, the shape graphs generated from the HCP data also contained: (1) a hybrid mesoscale structures of modularity as well as core/periphery; (2) a similar pattern of significantly less coreness scores for fixation (or resting) nodes as compared to cognitively demanding working memory nodes; and (3) temporal transitions extracted from TCM showing “within-task” transitions from one condition to the next (Supplementary Fig. 7).

As a third measure towards estimating the reliability of our approach, we tested the effect of parameter perturbation on shape properties and their relation to behavioral task performance. We varied both the TDA parameters—i.e., number of bins (or resolution, R) and percentage of overlap between bins (or gain, G)—to generate 49 different variations. Results are shown in Supplementary Fig. 8. Overall, the shape properties (e.g., the core-periphery structure) were reliably observed in a majority of parameter variations. Further, the association between modularity and task performance (in the CMP data) was also reliably significant across a majority of parameter variations. Altogether indicating robustness of results in the face of parameter perturbation.

Validation of the shape graphs against null models

To benchmark and validate the results generated from our approach, we employed three different null models. The first null model was designed to test whether the metrics generated from our approach (including the shape graph) were mainly driven by non-brain physiological signals (e.g., cardiac and respiratory) and spatiotemporal noise. In addition to these primary sources of noise, the individual variations in neuroanatomy were also included in the model, by using each anatomical scan (T1-weighted scan) for defining the baseline. The second and third null models were designed to test whether the task-unrelated non-stationarity primarily drove the shape graph and its properties. To generate last two null models, phase randomization (with and without constant phase sequence) of the original dataset was performed independently for each individual. After applying our approach to the data generated from these null models, we did not observe: (1) core−periphery or modularity characteristics in any of the null models, (2) task-dependent variation in the coreness scores, and (3) any structure in the temporal transitions across tasks (Supplementary Fig. 9).

In addition to these null models, we also tested whether head movement artifacts were responsible for individual differences in the TDA-generated shape graphs. No correlation was observed between the shape-graph properties (e.g., coreness score) and measures of head movement artifacts (e.g., mean relative and max absolute head displacement; all ps > 0.05).