Topological data analysis for revealing dynamic brain reconfiguration in MEG data

Ethics statement

This article utilised data collected for the HCP (Van Essen et al., 2012). The scanning protocol, participant recruitment procedures, and informed written consent forms, including consent to share deidentified data, were approved by the Washington University institutional review board (Van Essen et al., 2012). The IRB approval number is FEB-20220810-13614.

Subjects and data

The data we use in this research is the human non-invasive resting state and task Magnetoencephalography (MEG) data set which is publicly available from the Human Connectome Project (HCP) consortium (Larson-Prior et al., 2013). It is acquired on a Magnes 3600 MEG (4D NeuroImaging, San Diego, CA, USA) with 248 magnetometers and 23 reference channels at the sampling rate of 2034.5101 Hz. The data is available for a total of 100 subjects each performing three experimental paradigms: Sensory-motor, Working memory, and Story/Math.

The tasks in the sensory-motor paradigm involve the execution of a simple hand or foot movement. Which limb on which side is instructed by a visual cue, which serves to pace the movement. The working memory paradigm is similar or identical to the corresponding task acquired during fMRI imaging of HCP. Here, the participants have to remember the occurrence of n-back previously shown item (with n = 0 and n = 2) with the items being either tools or faces. Data are segmented to the onset of the non-target item (WM task). The story/math paradigm is the same as that is used in the fMRI component of the HCP (Barch et al., 2013). Participants listen either to auditory narratives of around 30 s duration or matched-duration simple arithmetic problems followed by a 2-alternative forced choice question (Binder et al., 2011). Subjects respond by right-hand button press (index or middle finger).

For every subject, all the paradigms consist of two experimental runs. A run consists of blocks of tasks, a block consists of several trials of the same task with a fixation period between the trials, and finally, a trial consists of a baseline and a stimulus. We shall also note that not all data is available for each subject.

For a single subject, the data acquired by Magnes 3600 MEG is processed as follows:

• Noisy channels with a high variance ratio and correlation to neighboring channels are detected and removed from further analysis.

• The bad channels and segments are removed with iterative independent component analysis (ICA) using spatial and temporal criteria (Mantini et al., 2011).

• Using ICA, independent components (ICs) are classified as ‘brain’ or ‘noise’ using six parameters: correlation between IC signals, the correlation between power time courses, the correlation between spectra, and three additional parameters derived from both spectral and temporal properties. Physiological artifacts are identified as magneto- and electro-cardiogram, eye movements, power supply bursting, and 1/f-like environmental noise. The details of this step can be found in Larson-Prior et al. (2013) and Mantini et al. (2011).

We note here two of its stages that are related to our preprocessing explained in the next section. First, the sampling rate is lowered to 506.6275 Hz, and second, the data from noisy channels are removed. As a result, the time points are reduced by 25% and the channels across the two runs might not be identical.

Preprocessing

We process the data further per subject and per paradigm. We lower the sampling rate down to 256 Hz. We remove all the time points corresponding to the fixation trials and the baseline keeping only the stimulus ones as well as the time points with missing values. Moreover, to concatenate the data from the two runs, we also remove the non-common channels across the runs. As a result, a subject for a given paradigm is represented by a matrix of the form (# channels) ×(# time points). The number of channels is changing between 200 to 248 depending on the subject, and the number of time points is ranging between 350,000 and 400,000.

At this point of the preprocessing, each subject is associated with three matrices, one for every experimental paradigm. As we want to compare the paradigms pairwise, in the next stage we concatenate these matrices two by two. This way, each subject is still associated with three matrices corresponding to the cases:

• Working memory and story/math;

• Working memory and sensory-motor;

• Story/math and sensory-motor.

These matrices are concatenated in the above order across the common channels in the paradigms. Even though the concatenation results in the loss of some channels, the loss is not significant and the resulting concatenated matrix has still 200 to 250 channels. In the meantime, the number of time points in the concatenated matrix is almost doubled and ranges between 700,000 and 800,000. In the end, the concatenated matrix is transposed so that the time points are on the vertical axis and the channels are on the horizontal axis.

Finally, we note that this vectorization of the MEG images causes the loss of the locations of channels relative to each other. However, it is shown in various studies that this process does not affect the success of the machine learning methods (Bray, Chang & Hoeft, 2009).

Topological data analysis: Mapper

Most of the network neuroscience studies utilize simple graphs focusing on dyadic connection ignoring higher-order interactions that could be crucial to extract insight across multiple scales (Torres et al., 2021). The existence, quantification, and comparison of these higher-order interactions necessitate the use of more advanced mathematical structures that can be studied by topological data analysis(TDA), which uses techniques from algebraic topology and computer science to analyze data sets (Centeno et al., 2022). Analysis of these non-dyadic relations makes it possible to deal with the open problems in network neuroscience (Andjelković, Tadić & Melnik, 2020; Billings et al., 2021; Guo et al., 2021; Helm, Blevins & Bassett, 2021; Patania et al., 2019; Santos et al., 2019; Saggar et al., 2018).

In our study, we adopt a TDA method called Mapper, which is a successful structure discovery and visualization technique for the exploration of high-dimensional data. The resulting mathematical graph from Mapper is a highly compact representation of the complex data revealing insightful coordinate-free visualization. Introduced in 2007 (Singh, Memoli & Carlsson, 2007), its application to biological data sets includes but is not limited to disease association, RNA folding, viral evolution, and immunology (Chan, Carlsson & Rabadan, 2013; Nielson et al., 2015; Li et al., 2015). Hence, its application to brain dynamics is promising as shown in the earlier studies (Geniesse, Chowdhury & Saggar, 2022; Geniesse et al., 2019; Saggar et al., 2018; Patania et al., 2019).

The construction of a Mapper graph from a point cloud is illustrated in Fig. 1: (i) The first step is the choice of a filter which assigns one (or more) values to each data point in the point cloud. The filter values can be height, coordinate values, a measure of centrality, or output of any data mining algorithm such as PCA, SVD, SNE, etc. (ii) The next step is to cover all possible filter values with overlapping intervals (or regions depending on the dimension of the filter). Three color-coded intervals covering the range of the height function are shown in Fig. 1. (iii) Next, the points whose filter values fall in the same interval (or region) clustered using a clustering algorithm such as hierarchical clustering, k-nearest neighbor (KNN), or single linkage clustering. For clustering, one can use any metric including correlation, Euclidean, L₁ or L_∞ metrics. (iv) The nodes representing the clusters are finally connected by an edge if the underlying clusters have a non-empty intersection.

The parameters of the algorithm are the number of intervals (or regions), the overlapping percentage, the distance metric, and the clustering algorithm. A high number of intervals increases the number of nodes in the final visualization; hence, defeats the purpose of having a highly compact representation. While an increase in the overlapping percentage results in a high number of edges and increases complexity, a low overlapping percentage produces disconnected clusters and misses the information about variation in the data due to underlying continuous filter values.

We apply the Mapper algorithm to the MEG data from HCP. We use the open source KeplerMapper Python package (Veen et al., 2019) to generate Mapper graphs from the minimally processed data. Our goal is to trace the brain activation patterns of each participant during working memory, story/math, and motor tasks. As explained above, the data is concatenated pairwise before entering the Mapper algorithm. It is not uncommon to use concatenated data to estimate task-state functional connectivity and brain networks (Richiardi et al., 2011; Hsu et al., 2014; Freeman, Donner & Heeger, 2011; Liu et al., 2014; Mokhtari & Hossein-Zadeh, 2013; Zhu et al., 2017). Moreover, comparative studies show that functional connectivity for the concatenated data was both qualitatively and quantitatively similar to that of continuous data during rest (Fair et al., 2007; Gavrilescu et al., 2008; Cheng et al., 2015) and task (Zhu et al., 2017).

The input data to Mapper is the (# time points) ×(# channels) dimensional matrix prepared by the preprocessing explained in Preprocessing.

We choose the Euclidean metric as a similarity measure between the vectors in the time-space (i.e., column vectors) of the input matrix. The euclidean metric is suitable in this setting as the ranges and means of the data columns do not vary significantly (Nielson et al., 2015). The choice of Euclidean metric is also proven to be successful in the previous applications of neuroimage data: In Romano et al. (2014), Mapper with Euclidean metric applied on sMRI data reveals high and low functioning neuro-phenotypes within Fragile X Syndrome. In Kyeong, Kim & Kim (2017), the method on fMRI data identified two unique subgroups of ADHD using the Euclidean metric. In another fMRI study, a data-driven search for different metrics indicates that the Euclidean metric best localizes outcome measures (Madan et al., 2017).

In the next step, the similarity information determined by the Euclidean metric is transformed into a low dimensional representation using a non-linear filter called t-SNE (Hinton & Roweis, 2002), which maintains the local geometry existing in the original time-space unlike more conventional linear filters such as PCA (van der Maaten & Hinton, 2008; Saggar et al., 2018). Multivariate and non-linear characteristics of inter-regional interactions suggest the use of non-linear methods such as t-SNE (Reinen et al., 2018; DiCarlo, Zoccolan & Rust, 2012). The data which is first reduced into two dimensions by t-SNE is then divided into overlapping bins. Following the earlier practice (Lum et al., 2013), the time points in each bin are further clustered using single linkage clustering with the Euclidean metric, which is computationally more efficient compared to the other clustering methods, and which does not require an initial number of clusters.

The common practice in Mapper applications is to test a large grid of parameters (i.e., overlapping percentage and number of intervals) to find the most stable graphs. Even though the stability of the Mapper algorithm under various parameters was studied before under certain conditions (Carrière & Oudot, 2018; Kalyanaraman, Kamruzzaman & Krishnamoorthy, 2017), we analyze several parameters of the algorithm to ensure the reliability of our result.

Centrality and community structure analysis

The mathematical graphs obtained from Mapper can be investigated by focusing its structures on different scales. It is established through many applications that the intermediate (mesoscale) structures identify certain characteristics that the local scale analysis of nodes (or edges) and global level of summary statistics are unsuccessful to detect. In this article, we concentrate on centrality and modularity mesoscale properties.

Centrality analysis of a network reveals the most important nodes (or edges) based on a quantification of node-node or node-edge relationships. The centrality of a node can be perceived as the communication ability with the other nodes or the closeness to the other nodes (Estrada, 2011). The centrality structure provides a perspective to comprehend how the brain states evolve during the ongoing task. The nodes with high centrality in the Mapper graph contain the most common brain activation patterns during the task. Here, we utilized four centrality measures: degree, eigenvalue, betweenness, and closeness. Based on the visual evidence from the box plots Figs. 2, 3, and 4 and the similar findings in Saggar et al. (2018), we claim for any one of those four centrality scores that

1. The mean of the centrality score of the nodes dominated by the working memory paradigm time points is greater than the mean of the centrality score of the nodes dominated by the story/math paradigm time points.

2. The mean of the centrality score of the nodes dominated by the sensory-motor paradigm time points is greater than the mean of the centrality score of the nodes dominated by the working memory paradigm time points.

3. The mean of the centrality score of the nodes dominated by the sensory-motor paradigm time points is greater than the mean of the centrality score of the nodes dominated by the story/math paradigm time points.

Modularity is one of the most commonly used metric to detect and characterize the community structure of the networks. Detecting communities in brain networks, are useful to identify the sub-networks that correspond to specialized functional components (Sporns & Betzel, 2016). In this article, we use modularity as defined in Newman (2006).

Analysis pipeline

We summarize all the steps explained so far in the previous sections in Fig. 5.

Earlier studies (Saggar et al., 2018; Saggar & Uddin, 2019; Duman, Tatar & Pirim, 2019) show that the topological properties of Mapper graphs are robust by construction to parameter perturbations. To ensure the reliability of the statistical results, we tested our null hypothesis using 48 different sets of Mapper parameters. The domains for # intervals, overlap %, and # clusters are {10, 15, 20}, {30, 40, 50, 60}, and {5, 10, 15, 20}, respectively. The parameter space is chosen in a way that the resulting Mapper graphs are connected to make modularity and centrality calculation possible. Statistically significant and reliable results are obtained for the large portion of parameter values (see Results for details).

To have a better understanding, we also track the steps of the proposed analysis with the story/math and sensory-motor paradigm time points of the subject, 106521.

After concatenating the paradigms across the common channels, the subject, 106521, is represented by a matrix of dimensions (357, 371 × 232). The first 177,162 time points of the 357,371 time points belong to the story/math paradigm and the rest to the motor-sensory. The Mapper algorithm with the projected data by SNE and the parameters 10, 50, and 10 corresponding to the number of intervals, percentage of overlap, and the number of clusters in the single linkage clustering, respectively, outputs the Mapper graph in Fig. 6. It has 960 nodes and 13,889 edges. The nodes are either in black if the majority of the time points in that node belong to the story/math paradigm or in green color if the majority of the time points in that node belong to the sensory-motor paradigm. In 70% of the nodes the story/math paradigm time points are in majority and these nodes are on the periphery of the Mapper graph. The remaining 30% of the nodes are dominated by the sensory-motor paradigm time points which are placed mostly in the center of the graph.

Code accessibility

The code described in the article is available from the corresponding author upon reasonable request.

Centrality of the Mapper graphs

For every pairwise scenario, working memory vs. story/math, working memory vs. sensory-motor, and story/math vs. sensory-motor and for every subject, we calculate four different centrality scores, degree, eigenvector, betweenness, and closeness centralities of the nodes in the Mapper graphs. Moreover, to show that the results do not depend on the Mapper graph, we repeat the same experiment with 48 different sets of Mapper parameters. The Mapper parameter space is made of triples of the form (# intervals, overlap %, # clusters). The domains for # intervals, overlap %, and # clusters are {10, 15, 20}, {30, 40, 50, 60}, and {5, 10, 15, 20}, respectively.

We look at the pairwise results. In the working memory and story/math paradigms, there are 60 common subjects. Figure 2 shows that, regardless of the Mapper parameters, all four centrality scores clearly distinguish the nodes dominated by the working memory time points from the nodes dominated by the story/math time points. We want to confirm this visual difference through statistical tests. The most relevant statistical test to use is the paired t-test that compares the mean of the centrality scores of the nodes dominated by the working memory time-points (μ_W) and the mean of the centrality scores of the nodes dominated by the story/math paradigm’s time-points (μ_S). However, the paired t-test only works under the assumption that sample points follow a normal distribution. Therefore, we shall check numerically using the Shapiro–Wilk normality test or visually using probability plots that the differences in centrality scores of different paradigms satisfy this fact.

There might be Mapper plots with the differences in the centrality scores not following a normal distribution. For example, the distribution of closeness centrality scores of the nodes of the Mapper graph with parameters 10-40-5 is non normal distribution with p = 0.02 Shapiro Wilk normality test (see Fig. 7 for the probability plot). For such cases, we use the paired permutation test (Good, 2013) with the test statistic being the mean of the differences between the centrality scores. We shall clarify how the p value of the test statistics is calculated. First, we calculate the observed mean differences μ_obs between the working memory and the story/math paradigms. The permutation test assumes that there is no difference between the test statistics. The translation of this null hypothesis to our case is that there is no difference between the paradigms. Hence, we can create new data sets by swapping the centrality measures of the paradigms for the same subject. We compute the mean of the differences between the paradigms for every possible new data set and compare them to μ_obs. If the alternative hypothesis is that the means are not equal, then the p value is 2min(p_l, p_g) where p_g (resp. p_l) is the probability of the mean of the differences being greater (resp. less) than μ_obs. If the alternative hypothesis is that one of the means is greater (resp. less) than the other mean, then the p is the probability of the mean of the differences being greater (resp. less) than μ_obs. As in the other hypothesis testing, if the p value is smaller than the significance level of 0.05, we say that the null hypothesis is not likely to be observed and favor the alternative. In Table 1, we illustrate the permutation test on the closeness centrality scores of the Mapper graph with parameters 10-40-5. We show how the paired permutation test creates a permuted data set. Later, the test statistics which is the mean of the differences in the centrality scores of the paradigms working memory and story/math are calculated. Going back to hypothesis testing, to verify our claims, we start with the null hypothesis:

1. The mean of the centrality scores of the working memory nodes is equal to the mean of the centrality scores of the story/math nodes.

Based on the Figs. 2, 3 and 4, we expect H₀ to be rejected for all the Mapper parameters and for all centrality score types. Then, we verify which mean is greater by simply comparing the sample means. We look at the pairwise combinations one by one. To our expectations, in working memory and story/math combination, the p values in the columns μ_W = μ_S of the Table 2 being less than 0.05 indicate that we shall reject the null hypothesis in favor of its alternative which is that μ_W and μ_S are significantly different. Moreover, as μ_S < μ_W and halves of all the p values in Table 2 are also less than 0.05, we can deduce that the mean of the centrality scores of the working memory nodes μ_W is significantly greater than the mean of the centrality scores of the story/math nodes μ_S.

In the case of story/math and sensory-motor combination with 43 common subjects, by going through similar steps as above, we come to the conclusion using the test results in Table 3, that the mean of the centrality scores of the sensory-motor nodes μ_M is significantly greater than the mean of the centrality scores of the story/math nodes μ_S.

For the final combination of working memory and sensory-motor with 21 common subjects, the Table 4 shows that the p values of either the paired t-test or the permutation test under the null hypothesis are all greater than the significance level of 0.05 indicating that there is not enough evidence to reject the null hypothesis. Hence, we conclude for all types of centrality scores that mean the mean of the centrality scores of the working memory μ_W and sensory-motor μ_M are not significantly different.

The above discussions show that, when compared pairwise, the centrality score of the nodes dominated by the time points of the story/math paradigm is the smallest. We summarize our findings in the Table 5.

Community structure and response time

Subjects’ performances during the working memory and story/math experimental paradigms are scored and timed. In this section, we investigate if the modularity structure of the Mapper graphs is related to these measurements.

The distribution of modularity scores of all subjects are shown by box plots for different set of mapper parameters in Fig. 8. x-axis of the figure represents different set of mapper parameter. The interval parameter given by the first two numbers are either 10, 15 or 20. The following two numbers for the cluster parameter are either 30, 40, 50 or 60. The remaining numbers represents the overlap percentage which is either 5, 10, 15 or 20. Once we give a closer look at Fig. 8, we observe that the median of the modularity scores tends to

1. increase with the increase in the interval parameter as an increase in that parameter results in more nodes, hence more edges in the Mapper graph (see Fig. 9),

2. decrease with the increase in cluster parameter as an increase in that parameter results in a decrease in the number of edges (see Fig. 10),

3. non-decrease with the increase in the overlap percentage parameter as an increase in that parameter up to certain level results in a high number of edges and beyond that level less number of edges since some of the small nodes is absorbed by the larger ones. According to Fig. 11, that level is between 50% and 60% for our data set.

Moreover, the experiment results show that the modularity score of the Mapper graphs is negatively correlated to the reaction time. According to Table 6, the negative correlation is observed in all the Mapper graphs with different parameters. All correlation scores, with the strongest recorded at the parameters (10 − 40 − 5) (see Fig. 12), are weak and non-significant.

We also note that the strongest correlation scores are associated with the lowest clustering parameters as an increase in this parameter decreases the number of edges within the same interval which more likely contains similar tasks (see Fig. 13).

Other methods

We also visualize the time points before the Mapper algorithm. In Fig. 14, we see that the story/math paradigm time points overlap with the sensory-motor paradigm time points which indicates that a distance-based clustering would not distinguish them.