Multiscale Functional Connectivity: Exploring the brain functional connectivity at different timescales.

Abstract

Human brains exhibit highly organized multiscale neurophysiological dynamics. Understanding those dynamic changes and the neuronal networks involved is critical for understanding how the brain functions in health and disease. Functional Magnetic Resonance Imaging (fMRI) is a prevalent neuroimaging technique for studying these complex interactions. However, analyzing fMRI data poses several challenges. Furthermore, most approaches for analyzing Functional Connectivity (FC) still rely on preprocessing or conventional methods, often built upon oversimplified assumptions. On top of that, those approaches often ignore frequency-related information despite evidence showing that fMRI data contain rich information that spans multiple timescales. This study introduces a novel methodology, Multiscale Functional Connectivity (MFC), to analyze fMRI data by decomposing the fMRI into their intrinsic modes, allowing us to separate the neurophysiological activation patterns at multiple timescales while separating them from other interfering components. Additionally, the proposed approach accounts for the natural nonlinear and nonstationary nature of fMRI and the particularities of each individual in a data-driven way. We evaluated the performance of our proposed methodology using three fMRI experiments. Our results demonstrate that our novel approach effectively separates the fMRI data into different timescales while identifying highly reliable functional connectivity patterns across individuals. In addition, we further extended our knowledge of how the FC for these three experiments spans among different timescales.

[name=Morante, color=BrickRed]M [name=Morante, color=Magenta]m [name=Auth02, color=ProcessBlue]A2

[inst1]organization=Department of Electrical and Computer Engineering, addressline=Aarhus University, postcode=8200, country=Denmark

1 Introduction

Even at rest, the brain exhibits intricate yet highly organized neurophysiological interactions that control and regulate our body’s main cognitive and physiological functions (Honey et al., 2007; Fair et al., 2009; Poldrack et al., 2011). Furthermore, although the constant interaction between different neuronal networks -along with endogenous somatic interactions- dominates the evolution of brain dynamics, empirical evidence has shown that environmental factors also affect this process by modulating or inducing changes in neuronal activity (Fox et al., 2005; Friston, 2011). As a result, brain activity displays highly complex spatiotemporal dynamic (Bolton et al., 2020b; Lurie et al., 2020).

The complexity of brain activity can also be seen in the ability of the brain to integrate and process a vast amount of information at different timescales, as discussed by Preti et al. (2017). Additionally, neuronal communities often exhibit complex interaction regardless of their anatomical proximity (Honey et al., 2007; Xu et al., 2016). All these natural features render complex spatial and temporal dynamics, which lead to the emergence of sophisticated cognitive processes such as perception, attention, or memory (Carpenter et al., 2021) and, ultimately, shape our thoughts, behaviors, and experiences, even the way we perceive and interpret our past (Schacter and Coyle, 1995).

Researchers have used a variety of neuroimaging techniques to study brain activity. Among all those techniques, functional Magnetic Resonance Imaging (fMRI) stands as one of the most popular for both neurophysiological, e.g., (Ozcelik and VanRullen, 2023; Morante et al., 2021; Chen et al., 2020; Chatzichristos et al., 2020), and clinical research, e.g., (Peng et al., 2021; Hannanu et al., 2020; Seo et al., 2019). Unlike other alternative noninvasive neuroimaging techniques, such as Electroencephalography (EEG), e.g., (Singh and Krishnan, 2023), or Magnetoencephalography (MEG) (Koshev et al., 2021), fMRI records brain activity indirectly by measuring the variation in the oxygenation levels on small volumes of tissue, referred to as voxels (Power et al., 2011), due to the metabolic oxygen consumption of the neurons, which is a process referred to as Bold Oxygenation Level-Dependent (BOLD) contrast Poldrack et al. (2011). Furthermore, the fMRI data contains a mixture of several interfering signals beyond brain-induced activity, such as movement, respiratory and cardiac pulsations (Bianciardi et al., 2009; Logothetis and Wandell, 2004).

The majority of fMRI studies, especially in the early days of fMRI research, have focused on identifying both the spatial areas associated with brain activity and their corresponding activation patterns, often in response to specific stimuli or tasks, usually referred to as sources (Friston, 2011; Fair et al., 2009). Methods such as the General Linear Model (GLM) (Poldrack et al., 2011), and matrix factorization techniques, such as Independent Component Analysis (ICA) or Dictionary Learning, stand as some of the most relevant conventional techniques (Morante, 2021). Nonetheless, those simple models are limited, as they lack the capability of exploring more complex interactions among neuronal networks and their Functional Connectivity (FC) (Bolton et al., 2020b).

In consequence, during the last two decades, the study of the FC and dynamic interaction of neuronal networks using fMRI has steadily gained more attention (Lurie et al., 2020), as Friston (2011) envisioned back in 2011. Researchers have steadily proposed a wide range of techniques to reveal the intricate nature of the brain’s FC (Bolton et al., 2020b; Lurie et al., 2020). Many of the most recent approaches – greatly inspired by the success of Machine Learning (ML) and Artificial Intelligence (AI) in a wide range of complicated problems – leverage machine or deep learning-based models due to their potential to uncover hidden interactions among brain networks that are not directly observable from the data (Kim et al., 2021).

No doubt, ML and DL can provide novel and powerful tools that appear well-suited for capturing FC and brain behavior. The literature (Bolton et al., 2020b; Pervaiz et al., 2020), however, shows that most of these approaches still use conventional preprocessing steps or rely on classical fMRI analysis techniques, such as ICA, as a critical preliminary step in their application. Therefore, even if all those novel cutting-edge techniques have great potential, their application will be hindered by the inherent assumptions and limitations of these preliminary steps.

For instance, conventional motion correction approaches rely on the parametric rigid-body model, which assumes that the brain only undergoes rigid motion. These approaches provide a reasonable estimation for main motion effects. Nonetheless, these motion components are richer than this simple model, as they can carry more information regarding certain behavioral and physiological aspects of the participants, such as arousal(Gu et al., 2019), and even age, as shown by Bolton et al. (2020a). Furthermore, these simple models overlook other more subtle sources of motion, such as respiratory or cardiac pulsations. These additional motion components are often challenging to identify since they appear highly structured, individual dependent, and even mimic other brain sources, as evidenced by Chen et al. (2020).

Temporal filtering also constitutes another relevant preprocessing step. From a general perspective, temporal filtering aims to remove portions of the frequency spectrum from the acquired fMRI data without relevant neurophysiological information (Poldrack et al., 2011). However, filtering is not without problems as it has been discussed in the literature (Lurie et al., 2020), even considering it inappropriate (Yuen et al., 2019). For example, several authors have pointed out that low-pass filtering can be problematic since it may remove signals of interest (Bolton et al., 2020a), while others have pointed out that high-pass filtering also has other downsides Lurie et al. (2020). In addition, the choice of filter parameters, such as cutoff frequencies, which can impact the analysis outcomes, is unclear since the optimal parameters vary among individuals. Consequently, finding an adequate filtering process that balances preserving relevant information and removing noise remains challenging in fMRI research.

On the other hand, on top of the standard preprocessing steps, many approaches for FC analysis utilize some classical fMRI source separation techniques. One of the most relevant approaches is ICA, which plays a crucial role in many analysis pipelines. In addition to separating neuronal networks of interest, ICA is widely used for separating interfering components, e.g., cardiac or motion artifacts, from neurophysiologically relevant sources or as a data-driven way to spatially parcellate the brain into relevant brain sources. Formally, ICA is a statistical-based method that assumes that the fMRI data can be described as a linear combination of a set of maximally independent sources (Sergios Theodoridis, 2020). Additionally, although often overlooked (Lenton et al., 2023), ICA also assumes that the fMRI is stationary from a statistical point of view.

Nevertheless, although there exist alternative approaches that rely on similar assumptions, there is compelling evidence of the brain’s nonlinear and nonstationary nature. Guan et al. (2020) explicitly studied the stationary and nonstationary behavior of different brain areas, showing that the brain dynamics are nonstationary and nonlinear. Furthermore, Zalesky et al. (2014) and Preti et al. (2017) reported consistent evidence supporting brain nonstationary behavior. As emphasized by Guan et al. (2020), researchers should consider nonlinear and nonstationary effects when studying brain FC.

Consequently, it becomes evident that there is a need for more advanced and adaptive preprocessing methods capable of incorporating the nonlinearity and nonstationarity nature of fMRI data. Considering all unique challenges brain activity poses, as we have discussed (Bolton et al., 2020b; Guan et al., 2020; Cordes et al., 2001), we propose a novel approach that separates neuronal brain activity from other interfering components and, at the same time, unravels those brain activity across different timescales, while naturally accounting for the dynamic, nonlinear, and multifaceted nature of brain activity, in a fully data-driven way tailored to each individual.

Put succinctly, our proposed methodology, referred to as Multiscale Functional Connectivity, relies on Multivariate Mode Decomposition (MMD) a Signal Processing technique that allows us to capture brain activity patterns at different timescales and separate them into neurophysiological components and interfering signals. To achieve this, MMD assumes that the fMRI data constitute a mixture of multiple inherent oscillatory components. Then, we use the results from this decomposition to study FC across different timescales.

In this regard, we would like to mention that the idea of decomposing fMRI data into their inherent oscillations has been explored by Yuen et al. (2019). However, unlike Yuen et al. (2019), who performed their study voxel-wise and utilized a band-pass filter for eliminating interfering components, we take some steps further by (a) incorporating the multivariate nature of fMRI, (b) avoiding any band-pass filter to separate the relevant components, and (c) exploring the FC associated with these oscillatory components further.

In summary, we propose an alternative way for extracting FC information from fMRI data, bypassing the limitations of conventional preprocessing steps by computing the nonlinear, nonstationary, and multivariate nature of fMRI data. In addition to aligning with the expected natural complexities of the fMRI data, our proposed approach can detect and separate interfering components, such as structural noise from heart and respiration components, and simultaneously offers new insights into how the brain networks operate at different timescales. Finally, we evaluated the performance of our proposed methodology using three distinct fMRI studies. Our results demonstrate that our novel approach effectively separates the fMRI data into different timescales while identifying functional connectivity patterns with high reliability across individuals. Furthermore, we further extended our knowledge of how the FC for these three experiments spans among different timescales.

2 Materials and Methods

2.1 Natural frequencies in fMRI data

To fully appreciate the motivation behind the proposed study, we need a clear understanding of the frequency behavior of the fMRI data. Firstly, while studies focusing on frequency-related aspects of fMRI are relatively sparse, existing research offers valuable insight into the frequency organization of the fMRI signal and brain dynamics. For instance, Cordes et al. (2001) demonstrates that the frequency contribution to the correlation patterns spans several frequency bands. Similarly, Yuen et al. (2019) investigated the inherent frequency components across different brain locations –in a voxel-wise fashion– yielding similar findings.

Unlike other commonly used alternative neuroimaging techniques such as EEG or MEG, in fMRI, neuronal activity is indirectly measured through the Blood Oxygen Level Dependent (BOLD) effect (Power et al., 2011), which limits the observable neuronal activation dynamics of the brain to lower frequencies (Preti et al., 2017). Similarly, the relatively low sampling ratio utilized by fMRI scanners also imposes a challenge when studying the frequency contributions to fMRI signals (Power et al., 2011), as it restricts the maximum potential accessible frequencies and, at the same time, can introduce aliasing with physiological components (Cordes et al., 2001).

Despite these limitations, fMRI frequency components comprise a rich spectrum that covers several relevant frequency bands. First, very low-frequency oscillations, lower than 10 mHz (Power et al., 2011), correspond to trends, scanner instabilities, and motion residuals. Neurophysiological activation patterns resulting from neuronal activity appear within the range of 10 to 200 mHz  (Cordes et al., 2001; Yuen et al., 2019) emphasizing the significant contribution of this frequency band to fluctuations related to brain activity, which corresponds with the natural band dominated by the BOLD response.

Additionally, fundamental respiratory oscillations occur around 250 mHz, while the first harmonic of respiration appears around 500 mHz (Frank et al., 2001). Contributions from blood vessels and cerebrospinal fluid pulsations fall within 400 to 800 mHz band. Cordes et al. (2001) also noted that those high-frequency components exhibited significant structured correlations among different brain areas due to the distinct anatomical distribution of the cerebral blood vessels and ventricles. Similarly, they also pointed out that the cardiac pulsation can spread to lower frequencies due to aliasing, appearing as additional interfering structured components, which complies with the observations by Soon et al. (2021).

2.2 Multiscale Functional Connectivity

In this study, we propose a novel methodology for analyzing fMRI called Multiscale Functional Connectivity (MFC). The fundamental idea behind this approach is illustrated in Figure 1. Observe that the process consists of a few steps. First, we collected the fMRI data, and we used MMD to unveil IMs associated with that particular individual. Then, we used the multivariate IMs to unveil their corresponding FC patterns. The result is a set of FC patterns that extend across several frequencies, providing a multiscale representation of the FC.

Herein, we explore in detail the most relevant aspects of each step.

Figure 1: Summary of the proposed methodology.

Multivariate Mode Decomposition (MMD)

As we discussed in the introduction, one of the novel points of the proposed approach consists of MMD to extract the natural oscillation within the fMRI data. Formally, MMD is a Signal Processing model that assumes that a multivariate signal of interest accepts a representation as a linear combination of a set of a particular family of amplitude- and frequency-modulated (AM-FM) functions, with a well-defined instantaneous frequency at any given time instance among all the channels (Huang et al., 1998; Rehman and Aftab, 2019). In other words, these intrinsic oscillations, i.e., IM functions, behave similarly to harmonics that remain relatively close to a particular frequency, yet they are flexible enough to accommodate fluctuations in both amplitude and frequency in a data-driven way (Dragomiretskiy and Zosso, 2014).

Unlike other similar alternative studies, such as the approach proposed by (Yuen et al., 2019), MMD inherently exploits the multivariate nature of fMRI data, as Figure 1 illustrates. This is a critical advantage and differentiates it from other approaches. The multivariate nature of MMD allows each ROI to exhibit distinct functional behavior while sharing a common oscillatory behavior, i.e., a common central frequency. This flexibility accommodates a wide range of response changes. This behavior contrasts with conventional approaches, such as ICA, where each voxel/ROI activity is decomposed as a linear combination of a single common time activation.

However, performing MMD constitutes a challenging task, and several methods have been proposed to solve it. Among all the available approaches, in this study, we will focus only on two of the most popular: Multivariate Empirical Mode Decomposition (MEMD) (Rehman and Mandic, 2009) and Multivariate Variational Mode Decomposition (MVMD) (Rehman and Aftab, 2019). Nevertheless, although both methods aim to achieve the same signal decomposition in a fully data-driven manner, they are vastly different from an algorithmic perspective.

On the one hand, MEMD aims to obtain the intrinsic modes using a greedy iterative process. This process unfolds by iteratively averaging the maxima and minima envelopes until we obtain an IM. Although this process is relatively simple for univariate signals (Huang et al., 1998), determining local extrema from multivariate signals is challenging. Rehman and Mandic (2009) solved this problem by projecting the multivariate signal in different directions on hyperspheres. Then, the extrema of the projections are interpolated component-wise to provide the multidimensional envelopes of the signal, which are then averaged to yield the multivariate mean signal. Finally, this average is subtracted from the original signal as an IM function. This process continues iteratively until the residual is sufficiently small.

Observe that this algorithmic procedure causes MEMD to extract the highest frequency first, followed by lower successive frequency components. This characteristic makes MEMD particularly sensitive to noise. Furthermore, due to MEMD’s greedy nature, any mistakes introduced during the extraction of the high-frequency components will spread through the subsequent low-frequency components.

On the other hand, MVMD has the same aim as MEMD but offers better performance (Rehman and Aftab, 2019). Unlike other methods, MVMD aims to estimate all modes simultaneously. The objective is to find an estimate of narrowband functions from the data, where the center frequency of each mode aligns across all the multivariate components of the signal. This renders MVMD more robust to noise and provides a better estimate of the modes.

Functional Connectivity extraction

As illustrated in Figure 1, we can use the obtained multivariate IMs associated with each particular frequency band to uncover the FC at different time scales, providing a complete multiscale FC representation of the fMRI data. Nowadays, we can find a wide range of approaches for estimating the brain’s FC (Bolton et al., 2020b). For simplicity, we propose to use conventional Pearson’s correlation coefficient analysis to calculate the correlation between the time activation patterns within each IM. Note that we can do this because, unlike Yuen et al. (2019) who performed a voxel-wise study, we unveiled a multivariate representation of each IM. Note that this approach differs from static FC analysis, where the raw FC signals are analyzed directly from the ROIs, often after applying static temporal filtering, or from conventional dFC, such as sliding-windows approaches (Lurie et al., 2020).

2.3 Experiments description and fMRI data

In this study, we considered three experiments from the WU-Minn Human Connectome Project (Van Essen et al., 2013). Specifically, we selected the resting state, motor, and gambling experiments from the HCP repository. In each experiment, we randomly selected 100 healthy participants aged 22 to 35 years.

The first experiment was resting-state, where participants were instructed to remain as still as possible during the scan. We chose this experiment because resting-state data is widely used for FC analysis, often providing reliable results.

The motor experiment followed a standard block paradigm, where a visual cue asked the participants to tap their left or right fingers, squeeze their left or right toes, or move their tongue. Each movement block lasted 12 seconds and was preceded by a 3-second visual cue. Additionally, there were three extra fixation blocks of 15 seconds each. We chose this simple experiment because the activation patterns and neuronal networks involved are well studied (Yeo et al., 2011; Buckner et al., 2011; Turner et al., 2018; Morante et al., 2021), facilitating the evaluation of the results.

Last but not least, the gambling experiment followed a random block paradigm, where participants tried to guess if a randomly generated number between 1 and 9 was either higher or lower than 5. There are two main reasons why we investigated this additional task-related experiment. Firstly, like the motor task, this experiment has been well studied, making it easier to evaluate the results. The second reason is that, unlike the motor task, the gambling experiment’s paradigm is unpredictable, i.e., the guesses of the participants cannot be determined a priori. This randomness adds an unpredictable effect to the responses, increasing the variability in the data and allowing for a more robust and consistent analysis. Additionally, we expect the level of arousal and effort for this experiment to be higher than for the other two experiments, which may appear reflected in the performance of MMD or the obtained inherent FC patterns.

Preprocessing and regions of interest

We obtained the fMRI data directly from the HCP repository. The three used datasets were collected using a 3T scanner with a repetition time (TR) of 720 ms. The specific descriptions of the experimental procedures and acquisition parameters are detailed in the HCP imaging protocols. Finally, on top of the standard preprocessing pipeline already applied by the HCP (Barch et al., 2013; Van Essen et al., 2013), we further smoothed each brain volume with a 4-mm FWHM Gaussian kernel.

For this study, we divided the brain into several regions of interest (ROIs) using the Automated Anatomical Labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002). Although the AAL atlas maps the entire brain, we only analyzed cerebral regions, which resulted in 90 ROIs.

ROI Region Module 7/8 Frontal middle FP Fronto-Parietal 9/10 Frontal middle orbital 11/12 Frontal inferior opercular 13/14 Frontal inferior triangular 15/16 Frontal inferior orbital 59/60 Parietal superior 61/62 Parietal inferior 1/2 Precentral TP Temporo-Parietal 17/18 Rolandic operculum 19/20 Supplementary motor area 29/30 Insula 57/58 Post-central 63/64 Supramarginal 69/70 Paracentral lobule 79/80 Heschl 81/82 Temporal superior 83/84 Temporal pole superior 89/90 Temporal inferior 71/72 Caudate BG Basal Ganglia 73/74 Putamen 75/76 Pallidum 43/44 Calcarine Occ Occipital 45/46 Cuneus 47/48 Lingual 49/50 Occipital superior 51/52 Occipital middle 53/54 Occipital inferior 55/56 Fusiform 3/4 Frontal superior DMN Default Mode Network 23/24 Frontal superior medial 31/32 Cingulum anterior 33/34 Cingulum middle 35/36 Cingulum posterior 65/66 Angular 67/68 Precuneus 85/86 Temporal middle 5/6 Frontal superior orbital Lim Limbic 21/22 Olfactory 25/26 Frontal medial orbital 27/28 Rectus 37/38 Hippocampus 39/40 Para-hippocampus 41/42 Amygdala 87/88 Temporal pole middle 77/78 Thalamus Th Thalamus

Table 1: Summary of the selected ROIs organized according to the functional modules described by Parente and Colosimo (2020). The numeric label refers to the 90 ROIs associated with the AAL from (Tzourio-Mazoyer et al., 2002), where odd and even correspond to the right and left hemispheres.

Following the recent module-based network organization proposed by Parente and Colosimo (2020), we grouped these 90 ROIs into seven functional modules. Table 1 contains information regarding the ROIs selected from the AAL and its modules. For the selected 90 ROIs, we extracted the related time series using Nilearn toolbox. Finally, we removed the mean value from each ROI.

Additionally, for the motor task, we added an extra 5 ROIs for the analysis of the time courses associated with the different parts of the motor cortex. For extracting these areas, we used the same motor templates for separating these motor ROIs as, for example, the one implemented by Morante et al. (2021).

3 Results and discussion

We evaluated the proposed methodology in three different fMRI experiments. Firstly, we conducted a comprehensive study of two of the most popular MMD algorithms using fMRI data. Then, we evaluated the neurophysiological relevance of the IMs and their relation to neuronal activity. Finally, we analyzed the FC patterns associated with the FC, their reliability among participants. Additionally, we highlighted the further insights gained through this approach, emphasizing its differences, advantages, and limitations over existing methodologies.

3.1 Intrinsic Mode decomposition analysis and evaluation

The first step of our proposed approach (see Figure 1) consists of extracting the IMs from the fMRI data by performing MMD. In this study, we focused on the two most popular algorithms for this decomposition: MVMD (Rehman and Aftab, 2019) and MEMD (Rehman and Mandic, 2009). Overall, this section aims to illustrate the process of extracting the IMs, understand their physiological meaning relative contributions to the fMRI signal, and, at the same time, examine which algorithm provides better results when working with fMRI data.

MMD results using MVMD

Figure 2: Frequency and energy distribution associated with each mode using MVMD. The boxplot depicts the corresponding results among all the studied participants for the resting-state, motor, and gambling fMRI experiments. The colored area highlights the modes within the neurophysiologically relevant frequency band.

Figure 2 illustrates the frequency (a) and energy (b) distribution of the IMs using the MVMD algorithm for all participants. Firstly, we observed that the first IM exhibits a dominant frequency centered around zero. Similarly, this mode shows a higher relative energy contribution among all the participants. A closer examination of the activation patterns among different individuals revealed that this mode captures trends and low-frequency variations, including residuals and motion-interfering components. We also expected a high energy contribution of this component due to the nature of these low-frequency trends; these components spread over the whole brain, and, therefore, their contribution to the overall detected signal shall be higher.

In addition, although we observed a high degree of consistency in frequency and energy among participants in this mode, a close observation of the actual activation patterns associated with each mode revealed that those trends are very different among participants, which indicates that those low-frequency trends are individual-dependent.

Modes 2 to 5, highlighted within a shaded area in Figure 2, appeared within the neurophysiologically relevant frequency range. This range includes typical brain activity frequencies related to several cognitive and neurophysiological processes (Cordes et al., 2001). Table 2 shows the average central frequencies for each IM among participants, and we found that those results were relatively consistent among participants and experiments.

Resting state Motor Gambling Freq. BW (mHz) Freq. BW (mHz) Freq. BW (mHz) Mode 1 0 6.8$\pm$ 1.0 0 21.6$\pm$ 1.5 0 26.3$\pm$ 1.5 Mode 2 23 14.3$\pm$ 1.9 23 36.5$\pm$ 4.1 26 44.9$\pm$ 4.3 Mode 3 54 15.0$\pm$ 2.2 56 40.5$\pm$ 1.9 66 50.8$\pm$ 2.3 Mode 4 96 15.9$\pm$ 2.7 99 43.0$\pm$ 3.1 115 53.0$\pm$ 2.9 Mode 5 156 16.9$\pm$ 2.7 157 47.6$\pm$ 3.9 182 57.4$\pm$ 3.6 Mode 6 232 16.5$\pm$ 2.8 231 50.6$\pm$ 3.6 257 57.0$\pm$ 3.8 Mode 7 310 17.2$\pm$ 3.1 308 50.8$\pm$ 3.9 330 57.6$\pm$ 4.1 Mode 8 407 18.1$\pm$ 2.5 400 54.1$\pm$ 4.5 419 61.6$\pm$ 4.4 Mode 9 521 18.3$\pm$ 2.4 514 56.6$\pm$ 3.9 533 63.8$\pm$ 4.0 Mode 10 632 17.4$\pm$ 2.3 627 53.8$\pm$ 5.2 627 56.8$\pm$ 7.4

Table 2: Average frequency and their corresponding bandwidth (BW) associated with each IM for the different studied experiments for MVMD.

The remaining modes spread across high frequencies. As discussed in Section 2.1, signals within this frequency range originate from a mixture of different interfering components. These components include those induced by respiration movements, heartbeat, and cerebrospinal fluid pulsations. Notably, mode 6 –approximately centered at 250 mHz– captures the primary respiratory-related component, while mode 9, with a central frequency of approximately 520 mHz, aligns well with the first harmonic of the respiratory frequency (Cordes et al., 2001; Yuen et al., 2019).

Regarding the relative energy contribution, according to Figure 3 (b), which illustrates the relative energies of IMs among all participants. We observed that, in general, modes with lower frequencies exhibited higher energy, while those with higher frequencies had less and less energy contributions. Specifically, the first three modes contain the most energy, whereas the high-frequency modes contribute comparatively little.

Overall, we observed that the results obtained from MVMD align with the expected behavior of the data as described in the existing literature, as we have briefly summarized in Section 2.1. These results validate MVMD as a highly suitable method for this particular task, as it effectively captures the different IMs naturally present in fMRI among participants.

MMD results using MEMD

Figure 3: Frequency and energy distribution associated with each mode using MEMD. The boxplot depicts the corresponding results among all the studied participants for the resting-state, motor, and gambling fMRI experiments.

Figure 3 illustrates the frequency (a) and energy (b) distribution of the modes derived from analysis using MEMD among all the participants. However, in contrast to the results of MVMD shown in Figure 2, the results from MEMD exhibit an entirely different behavior.

First, we observed that the number of intrinsic modes recovered by MEMD varies among individuals and experiments. Figure 3 shows that the maximum number of IMs among experiments is also different per experiment. No doubt, variations in the number of IMs per participant and experiment are expected; it is natural to assume that not all individuals activate the exact number of IMs, particularly if they are performing different tasks. In this regard, MEMD offers more flexibility by using a greedy approach and, unlike MVMD, which fixes the number of IMs to a specific value, MEMD finds all potential components without the necessity of setting any prior guess. However, although this feature is desirable in practice, the greedy approach makes MEMD particularly sensitive to noise. In this way, the results obtained from MEMD are poor and seem rather unusual due to several factors we discuss below.

On the one hand, we observed that the resting-state experiment was the experiment that exhibited the maximum number of IMs, displaying at most 11 modes. In contrast, gambling and motor only exhibited a maximum of 10 and 9, respectively. This unexpected result contradicts the expected nature of these experiments; in the resting state, participants were at rest and calm. Consequently, we would have expected to find the minimum number of components when the participants rested. On the other hand, participants were more actively engaged during the performance of the gambling experiments. Therefore, we would have expected to find the maximum number of active IMs associated with this experiment instead.

Similarly, the overall frequency distribution is odd. We found that the modes appeared predominantly in the low-frequency range. For instance, all the modes, except the last two from the resting-state experiment, are below 100 mHz. Table 3 shows the central frequencies and bandwidths associated with these modes. From observing this table, it is evident that many modes appeared to be dominated by the low frequencies and many of those modes show overlapping bandwidths.

Resting state Motor Gambling Freq. BW (mHz) Freq. BW (mHz) Freq. BW (mHz) Mode 1 1 4.6$\pm$ 0.1 5 19.4$\pm$ 0.3 5 24.1$\pm$ 0.7 Mode 2 2 5.1$\pm$ 0.8 6 20.6$\pm$ 4.6 9 27.3$\pm$ 5.4 Mode 3 3 6.7$\pm$ 1.2 12 27.6$\pm$ 5.0 18 37.4$\pm$ 7.3 Mode 4 6 8.5$\pm$ 1.1 23 35.4$\pm$ 4.3 34 46.5$\pm$ 4.6 Mode 5 11 9.9$\pm$ 1.2 44 41.6$\pm$ 5.4 63 54.4$\pm$ 4.9 Mode 6 19 12.1$\pm$ 1.9 81 51.9$\pm$ 8.3 109 68.2$\pm$ 10.2 Mode 7 32 14.5$\pm$ 2.4 142 60.6$\pm$ 9.3 202 78.3$\pm$ 12.7 Mode 8 51 14.4$\pm$ 3.3 192 60.2$\pm$ 10.5 222 80.2$\pm$ 13.3 Mode 9 82 12.8$\pm$ 3.4 190 62.0$\pm$ 5.4 244 76.2$\pm$ 19.2 Mode 10 140 12.6$\pm$ 2.5 170 52.8$\pm$ 9.6 Mode 11 153 12.3$\pm$ 2.5

Table 3: Average frequency and their corresponding bandwidth (BW) associated with each IM for the different studied experiments for MEMD.

We were surprised by not observing any high-frequency component related to respiratory or cardiac pulsations. High-frequency contributions from physiological sources, including respiratory components, as well as from cerebrospinal fluid and cardiac pulsations –as we discussed in previous sections– constitute a natural part of the fMRI signal. Somehow, these contributions are absent. Overall those results contrast strongly with the results from MVMD in Figure 2, which was able to capture the full range of frequencies naturally associated with fMRI data.

Finally, when analyzing the energy distribution in Figure 3, we observed that the first mode in all the experiments exhibited the most energy contribution, as expected due to the nature of the first mode, which is in line with the discussion of the MVMD results. However, unlike the results in Figure 2, we observed that the energy associated with the rest of the IMs from MEMD exhibited almost a homogeneous distribution.

We hypothesized that MEMD led to those unnatural and unreliable results due to the combination of two factors. First, the nature of the noise in the data, fMRI data is corrupted with Rician noise Poldrack et al. (2011), with especially high levels of energy at low frequencies due to motion artifacts. Second, the greedy nature of MEMD means that the algorithm is mostly driven by the noise associated with the first retrieved components rather than the inherent oscillation of the data.

In addition, another potential problem that might have influenced the results of MEMD is the multivariate extension of this algorithm. As illustrated by Yuen et al. (2019), the univariate version of these algorithms, namely EMD and VMD, exhibited similar results. However, this is not the case we observed for the multivariate case. On the contrary, we observed that the results obtained from MEMD exhibit an entirely different behavior.

Overall, we observed that MEMD seems to be driven by low-frequency trends, as most of the components appear dominated by those low trends (see Table 3). Therefore, after all these observations, we conclude that MEMD appears unsuitable for MMD in fMRI data, as it is unable to extract reliable information from the data.

Frequency and energy differences between fMRI experiments

Based on the results from the previous section, it is clear that MVMD is an excellent candidate for MMD in fMRI data. Specifically, it provides a good separation of the intrinsic modes of the fMRI across the expected natural frequencies. This section focuses on the analysis of the results obtained from MVMD among the studied fMRI experiments.

Overall, the frequency and energy distribution for MVMD in Figure 2 show a similar trend among all the studied experiments. Nonetheless, upon closer examination, we noticed some interesting differences among experiments.

For instance, modes 5 and 6 in the gambling experiment exhibited a significantly higher central frequency than the other experiments. Similarly, when examining the energy distribution of these modes, we observe that modes 6 and 7 significantly contribute more to the overall energy of the signal.

The nature of mode 6 makes this result particularly interesting. These findings suggest that the participant’s physiological state during the gambling experiment may have been different from the resting-state and the motor task experiments. The higher frequency rate of the respiratory component associated with mode 6, as well as the high energy contributions in modes 6 and 7, which contain cardiac contributions, seem to indicate faster physiological rhythms. This particular physiological state corresponds with a higher level of arousal, which potentially reflects a state of excitement or nervousness from the gambling experiment, which contrasts with the potential physiological states from the other experiments.

Moreover, a closer examination of the energy distribution of low-frequency intrinsic modes among the different experiments reveals a significant difference in their energy contributions. In the resting-state experiment, the majority of energy is concentrated in the first mode, representing the dominant and most energetic component. This mode is followed by modes 2 and 3, which also exhibit significant energy levels. In contrast, in the task-related experiments, the first mode is less dominant, and neurophysiologically relevant modes (modes 2-5) consistently show higher energy levels. Modes 2 and 3 capture a substantial portion of the energy, while modes 4 and 5, although showing lower energy contributions, are higher compared to the resting state.

This observed energy distribution between IMs is consistent with the nature of the different explored experiments: during the resting-state experiment, the participants were instructed to remain at rest, with their eyes closed, and think of nothing. Therefore, the brain activity at resting should have a relatively low contribution to the overall signal’s energy. By contrast, during task-related experiments, the brain actively participates in a particular task, using more energy as a consequence.

The difference in energy contributions of neurophysiological IMs between resting-state and task-related experiments suggests that the brain is more engaged and active during task performance compared to a state of rest, which makes intuitive sense. The higher energy levels in neurophysiologically associated modes indicate increased neural activity and processing during the task, highlighting the cognitive demands and involvement required for task completion.

3.2 Further analysis of the nature of the IMs

As Figure 1 illustrates, the MMD decomposition offers an alternative way to decompose the fMRI data in inherent modes. As illustrated in Figure 1, each IM provides a multivariate response where each ROI exhibits its particular response. However, despite this flexibility, all of those time activation patterns are characterized by sharing a common oscillatory behavior (Rehman and Aftab, 2019). Therefore, having such flexibility, it is still unclear how those different time activation patterns behave among ROIs, and it is unclear if they still provide reliable and consistent results among participants.

Consequently, to further shed light on these questions, we focus our attention only on the neurophysiologically relevant modes, i.e., IMs 2, 3, 4, and 5. For this study, we decided to investigate only the time activations using the task-related experiments because, unlike the resting-state experiment, for the task-related experiments, we have access to experimental paradigms, which allows us to approximately define the actual hemodynamic response of the brain and use it for comparisons.

Examination of the activation patterns for the motor task experiment

For the motor task experiment, we selected the primary visual cortex (ROI 43 and 44) for the visual responses. In addition, for this study, we further divided the ROI from the motor cortex into five additional motor-related ROIs for the different motor areas related to the right/left hands, feet and tongue. Finally, we performed MMD again with these additional ROIs.

Figure 4 illustrates the average time courses for the relevant ROIs associated with the motor experiment among all the participants. The reported lines correspond to the average time activations for modes 2 and 3 among all the studied participants (without any additional postprocessing). Finally, the orange lines in Figure 4 represent the canonical task-related component expected within each main ROI. We obtained the task-related components using the classical convolutional model with the canonical Hemodynamic Response Function (HRF) (Power et al., 2011). This visualization is crucial because it reflects MVMD’s capability to reveal neurophysiological information, and how each mode contributes to the brain activity.

Figure 4: Average time courses among all the studied participants associated with the first neurophysiological modes (2 and 3) from the results of MVMD (blue lines), and the canonical task-related component (orange line) for the specified ROIs associated with the different motor tasks.

Those results showed that IM2 and IM3 effectively capture information related to the expected brain activation patterns within their corresponding ROIs of interest. For instance, mode 2 closely aligns with the block-related activity from the motor cortex ROIs. On the other hand, mode 3, which exhibits a higher central frequency (see Table 2), encodes the visual cue associated with the motor task and fixation.

Interestingly, although mode 3 seems better suited to visual cues, we can see that it also contributes to the motor cortex, which indicates a potential connection between the visual and motor-related areas. This is a very interesting feature because this allows us to we can explore how the different areas interact at different timescales.

Overall, these findings provide strong evidence that MMD is an effective method for analyzing brain activity and extracting meaningful information from different regions of interest. The fact that mode 2 aligns closely with block-related activity in the motor cortex ROIs, while mode 3 captures the visual cues associated with the motor task and fixation, demonstrates the interpretability and naturalness of the results obtained through MMD. This further supports the validity and usefulness of MMD in studying the connection between different brain regions and their functional roles.

3.3 Multiscale Functional Connectivity

Reliability of the FC maps among participants

Reliability refers to the ability of methods to consistently detect significant activity within the expected ROIs (Morante et al., 2020). In this case, we focus the attention on the reliability associated with the FC patterns from MMD among participants. Since all participants are analyzed independently, failing to provide consistent results among participants indicates that the method is unable to capture relevant common information.

Therefore, we conducted a study to assess the reliability of the FC among all participants. This study aimed to demonstrate how different FC patterns behaved across participants. In particular, we studied the individual reliability of the static FC patterns associated with each IM for the three considered fMRI experiments.

For completeness, we evaluated the FC across all IMs. To this end, we calculated the Pearson’s correlation obtained from all the possible pairs of comparisons across all participants. We want to emphasize that this step was critical in ensuring the validity and generalizability of our findings, as it allowed us to confirm that the FC patterns associated with each IM were consistent and reliable across all participants, providing a solid foundation for our subsequent analyses.

Figure 5 illustrates the similarity of the FC associated with each mode using MVMD for the three studied experiments. The boxplots depict the results obtained across all pair comparisons for all participants.

Figure 5: Individual reliability of the static FC patterns associated with each intrinsic mode for the three considered fMRI experiments using MVMD. The boxplots depict the Pearson’s correlation values obtained among all the possible individual pair comparisons across all the participants

Upon examination, we observed that all experiments followed a similar trend. Overall, modes within the neurophysiological activity band exhibited higher similarity than the other modes. Specifically, modes 3 and 4 demonstrate remarkable reliability.

In contrast, mode 1 exhibited the lowest similarity across participants in all experiments. These results indicate that the observed patterns associated with the first mode largely depend on each individual; as we already pointed out in the previous section, the first mode contains individual trends and motion residuals that depend on each particular individual.

Regarding the remaining modes above the neurophysiological band, even though they are less consistent than modes within the neurophysiological band, they still exhibit a higher similarity to mode 1. This result indicates that the FC associated with these IMs exhibits a certain degree of consistency among individuals, despite mostly containing interfering physiological components as previously discussed.

Although the three studied experiments exhibited a similar trend, a closer examination reveals some relevant differences. First, the resting-state experiment consistently showed increased similarity for all modes, with mode 4 (96 mHz) being the most consistent among all individuals. This may also be related to why resting-state is the dominant approach in functional connectivity (FC) studies; on top of its simplicity, these results highlight that, in contrast to task-related experiments, activation patterns at different frequencies remain stable and produce more reliable FC patterns among many individuals.

Interestingly, if we focus on the high-frequency modes, modes 8 and 9 maintain a relatively higher similarity than modes 7 or 6. We speculate that this is due to cardiac interfering components at rest, which provide similar patterns. Cardiovascular and cerebrospinal pulsations during rest are more stable and predictable, leading to consistent connectivity patterns across individuals in the high-frequency modes.

Conversely, for the motor task, we observed that mode 3 (56 mHz) was the most consistent among all the participants, followed closely by mode 4. This result was expected since mode 3 (56 mHz) contains a visual cue common to all participants. The presence of a shared stimulus in this mode likely contributes to its higher consistency among individuals during the motor task.

Last but not least, the gambling experiment displayed the lowest similarity across all modes. Nonetheless, we still observed that modes 3 and 4 exhibited the highest consistency among the other modes. Interestingly, we noticed a significant drop in the similarity in the modes at the higher frequencies.

For completeness, we further evaluated the reliability of the modes obtained using the MEMD algorithm. Figure 6 illustrates the obtained results using MEMD.

Figure 6: Individual reliability of the static FC patterns associated with each intrinsic mode for the three considered fMRI experiments using MEMD. The boxplots depict the Pearson’s correlation values obtained among all the possible individual pair-comparisons across all the participants

No doubt, we observed that MEMD results exhibited an entirely different behavior. In contrast to the results from Figure 5, we noted that the first modes of MEMD are poorly similar among participants, whereas the high-frequency modes exhibited the highest similarity. As we previously hypothesized, these results support the idea that the MEMD algorithm seems to be mostly driven by the noise associated with the first retrieved components rather than the inherent oscillation of the data.

Overall, we can see that those results are difficult to interpret and are inconsistent with the expected behavior of the data. Therefore, these results further evidence that MEMD is unsuitable.

MFC patterns and analysis

After understanding the frequency and energy distribution of the different IMs, as well as their reliability among individuals, we focus on the analysis of the FC patterns associated with the IMs. In this case, we are interested in studying only neurophysiological activation patterns. Therefore, we analyzed only the FC patterns associated with modes 2, 3, 4, and 5.

Figure 7 depicts the average FC patterns for each mode. We obtained those FC patterns by averaging the individual FC patterns across all participants. Each row corresponds to a particular mode, and each column contains different experiments. For all comparisons, we performed a statistical test with respect to a null dataset generated from each particular decomposition by randomly mixing the temporal samples of the IMs. Pearson’s correlation coefficients were Fisher-Z transformed. The lower diagonal of each connectivity matrix displays the average correlation coefficients, while the upper diagonal shows only the correlation values that also exhibited significant activation compared to the null data derived from a permutation-based t-test corrected to false discovery rate adjusted to $p<0.001$. For convenience, we arranged the ROIs according to the leading module (left and right), following the order reported in Table 1.

Figure 7: FC patterns for the modes 2, 3, 4 and 5 for the three studied fMRI experiments. The mean FC patterns were estimated by averaging across 100 participants. Person’s correlation coefficients were Fisher-Z transformed. The lower diagonal part shows all the averaged correlation coefficients. The upper diagonal only displays significant correlation coefficients compared to the null dataset from a permutation-based t-test corrected with a false positive rate adjusted to $p\leq 0.001$.

Overall, Figure 7 indicates similarities among the experiments, particularly within the networks linked to the main modules. A significant activation pattern was observed in modes 2, 3, 4, and 5. On closer inspection, we observed some relevant patterns. First, mode 2 was characterized by extensive activity across several brain networks involving numerous ROIs.

This connectivity suggests a broad engagement of brain regions in this mode, indicating a complex and integrative role in coordinating diverse neural processes. In contrast, Mode 5 demonstrated more focused activation patterns, indicating less extensive networks. This mode’s more localized activity suggests specialized functions, potentially related to specific cognitive tasks. Further visual inspection of the FC matrices showed that the strength and distribution of the correlations within and between brain modules varied significantly across the different modes and experiments.

For instance, modes 3 and 4 exhibited similar connectivity patterns. These modes contributed to overall functional connectivity by integrating and modulating activity across. In mode 2, a notable correlation was observed between the occipital and temporal modules. This connection, however, was considerably weaker in the motor and gambling experiments, suggesting contextual modulation of connectivity. Interestingly, while the temporal cortex and the Default Mode Network (DMN) were disconnected in the resting-state experiment, they were conspicuously linked during task-related activities. This finding highlights brain connectivity dynamics and adaptability across different tasks.

Furthermore, during the gambling experiment in Mode 2, a pronounced left-hemisphere lateralization was linked with the limbic module, showing the potential role of emotional processing in gambling tasks. In contrast, the DMN and the Occipital cortex exhibited a significant correlation during the motor task, which was absent in other experiments. This suggests a specific involvement of these regions in motor control and visual processing.

Mode 3 and 4 mirrored the connectivity patterns of mode 2 but with more intense connections between the occipital and temporal cortex, indicating a potential enhancement or amplification of sensory processing in this mode. Unlike modes 2 and 3, which did not show the same structural pattern, hinting at functional differences even within those modes.

Mode 5 presented a unique scenario; it exhibited almost identical patterns across all experiments but with noticeably reduced connections both within and between modules. Notably, significant activity in the limbic module was observed in its connections with other modules, except in the gambling experiment. This absence might reflect the unique emotional or cognitive demands of the gambling task, which could dissociate limbic activity from other brain regions (Harrison et al., 2021).

In modes 2 and 3, the consistent significant correlation between the occipital and temporal modules indicates strong sensory integration. This is less pronounced in motor and gambling experiments, suggesting a shift of resources towards task-specific networks. Occipital-temporal connectivity is critical for processes such as visual recognition and memory retrieval, which may be differentially recruited depending on the task at hand.

The temporal cortex and DMN’s distinct connectivity in task-related experiments, as compared to the resting state, underlines the task’s influence on the brain’s intrinsic activity. The DMN, typically active during rest and internal mentation, plays an interactive role when engaging with external tasks. This is reflected in enhanced temporal-DMN connectivity.

4 Limitations and future work

As with any other study, it is important to take into account there are some limitations concerning our study.

Firstly, our analysis was conducted on three fMRI datasets with similar imaging protocols. While these datasets provided valuable insights into MFC behavior and performance, it is important to recognize that the generalizability of our findings may be limited. To fully understand the robustness and applicability of MFC, it would be beneficial to apply this methodology to a broader range of datasets, conditions, and experiments.

Secondly, for simplicity, our study focused on the Region of Interest (ROI) level. This decision was made to simplify the analysis and provide clear, interpretable results. However, the MFC methodology is inherently versatile and can be applied directly at the voxel level as well. Exploring the voxel-level application of MFC could yield fine-grained insights into brain functional connectivity patterns.

Last but not least, in this study, we did not explore in detail the dynamic aspects of the MFC to facilitate the introduction and understanding of the proposed approach. While our study provided preliminary insights, we did not explore completely all the advantages that the temporal time activation patterns of the IMs have to offer. A more thorough examination of temporal activation patterns, connectivity dynamics, and how they relate to conventional dynamic FC approaches is necessary to fully understand MFC capabilities.

Regarding potential paths and future work, there are several promising directions for future research to build upon our current findings and address the identified limitations.

One potential direction is to investigate MFC performance using more advanced algorithms tailored specifically to fMRI data. In our current study, we employed the two most popular algorithms, MEMD and MVMD, with quite different success. We believe that the exploration and development or adaptation of a more sophisticated algorithm could potentially improve the reliability and accuracy of the reported results, leading to a better characterization of the FC.

Another relevant path for future research involves exploring group-level strategies that integrate information from multiple participants. Our current study focused on individual-level analysis for MFC, which exploits MMD’s natural flexibility for accommodating individual variations. Nonetheless, it is still unclear how group-level strategies could be potentially integrated into MFC to accommodate for common information and, potentially enhance the robustness and generalizability of our findings.

5 Conclusions

In this study, we proposed a novel methodology for extracting neurophysiological functional information from fMRI data across multiple timescales, referred to as Multiscale Functional Connectivity (MFC), by decomposing fMRI activity into distinct intrinsic modes using Multivariate Mode Decomposition (MMD). To the best of our knowledge, this is the first time such an analysis has been performed on fMRI data. Unlike previous similar studies, such as Yue et al. (2019), our proposed method explicitly exploits the multivariate nature of fMRI and accommodates interfering components without manually specified temporal filters. As such, our method provides a more comprehensive view of the underlying fMRI data, which is essential for understanding its dynamics and interactions.

The analysis of three different fMRI experiments revealed that the algorithm Multivariate Variational Mode Decomposition (MVMD) constitutes a suitable candidate for MMD in fMRI data. Our results showed that MVMD provides a meaningful representation consistent with previous research while offering reliable, functional connectivity patterns among participants. The study of functional connectivity patterns among different neurophysiological components sheds light on the intertwined roles of various brain networks at several timescales. We believe that further adoption of MFC would give more insight into fMRI data and enhance our understanding of brain behavior.

Funding

This work has been partly supported by Lundbeck Foundation Investigator network grant R426-2023-133.