Formally, as introduced by Rehman and Aftab (2019), given a multivariate signal x(t) containing data from C channels, i.e., x(t)=[x1(t),x2(t),…,xC(t)]T, the MVMD algorithm tries to decompose the observed data as a linear combination of K principal multivariate oscillations, u^(k)(t), as in Equation 1, by solving the following optimization problem:

where ŭc(k) stands for the Hilbert Transform of uc(k).

Observe that this particular formulation sets multiple linear constraints, each one corresponding to each particular channel. On the other hand, the main loss function works over the K intrinsic modes, where are minimized to remain close to a common central frequency, ω_k.

Solving the optimization task in Equation 2 posses several challenges due to the inherent complexity of the task and the constraints involved. Nonetheless, (Rehman and Aftab 2019) resorted on the classical divide-and-conquer approach, and employed the Alternating Direction Method of Multipliers (ADMM) (Boyd et al., 2011) to effectively integrate multiple constraints. Put succinctly, this approach effectively solves the described optimization task by splitting the problem in multiple easier-to-solve subproblems. The resulting problems are then solved iteratively, producing an estimate of the intrinsic modes. See the original work by Rehman and Aftab (2019) for the particular formulation details and optimization steps.

Although MVMD is one of the most popular algorithms for performing MMD, several alternative methods have been proposed. These include Multivariate Empirical Mode Decomposition (MEMD) (Rehman and Mandic, 2009), multivariate iterative filtering (Cicone and Pellegrino, 2022), and multivariate chirp mode decomposition (Chen Q. et al., 2020). Among these, MEMD stands out as a widely used alternative to MVMD due to some unique features. For instance, MEMD does not require specifying the number of intrinsic modes beforehand (see Supplementary material, Section B for further details), which makes it attractive in cases when it is not possible to find a good estimate for the number of modes.

In essence, MEMD is a data-driven approach that decomposes a multivariate signal into its intrinsic oscillatory components through an iterative greedy process. Similar to MVMD, MEMD aims to extract the same intrinsic modes u^(k)(t) as described in Equation 1. However, MEMD employs a distinct iterative sifting process. At each iteration, it estimates a local mean by averaging the maximum and minimum envelopes of the signal. This local mean is then subtracted from the analyzed signal to produce an IM.

Despite its effectiveness in decomposing multivariate signals, MEMD has a significant drawback: it is highly sensitive to noise. This sensitivity can lead to a wide range of problems, including mode mixing, where a single mode contains multiple frequency components (Eriksen and Rehman, 2023). In addition, MEMD is particularly sensitive to noise, as high level of noise may interfere with the detection of the where ŭc(k) stands for the Hilbert Transform of uc(k). For all these reasons, MEMD does not appear as a suitable algorithm for fMRI analysis.

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.

Overall, fMRI frequency components comprise a rich spectrum that covers several relevant frequency bands. First, very low-frequency oscillations, much lower than 10 mHz (Power et al., 2011), correspond to trends,or signal drifts. Unlike low-frequency components (Tong et al., 2019), trends and drifts have been consistently been attributed to a combination of physiological fluctuations (Huettel et al., 2009), head motion residuals (Maknojia et al., 2019), and scanner instabilities (Smith et al., 1999). Neurophysiological activation patterns resulting from neuronal activity appear within the range of 10 to 200 mHz (Cordes et al., 2001; Yuen et al., 2019; Lurie et al., 2020) 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 the 400 to 800 mHz band. Similarly, 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 (Cordes et al., 2001). Similarly, they also pointed out that cardiac pulsations can spread to lower frequencies due to aliasing, appearing as additional interfering structured components, which complies with observations by Soon et al. (2021).

In this study, we propose a new method for analyzing FC using MVMD. The Figure 1 illustrates the main steps of our proposed approach:

In this study, we investigated three experiments from the WU-Minn Human Connectome Project (HCP) (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, with eyes open with relaxed fixation on a projected bright cross-hair on a dark background and presented in a darkened room³. 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 s and was preceded by a 3-s visual cue. Additionally, there were two 15-s fixation blocks. We chose this simple experiment because the activation patterns and neuronal networks involved are well documented (Morante et al., 2020), 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, this experiment has been well studied, making the evaluation process easier. The second reason is that, unlike the motor experiment, 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.

We examined the intrinsic modes obtained using MVMD among all the different experiments. Figure 2 illustrates the central frequencies, ω_k, (a), and relative energy contribution (b) of the IMs for all the studied participants. We calculated the relative energy contribution associated to each mode as the ration between the energy of the k-th mode, say E_k, and the total energy of the signal, E_tot, i.e., E_k/E_tot. Each colored box depicts the results associated with each studied experiment. For convenience, the shadowed area highlights the IMs within the neurophysiological bandwidth.

Regarding the intrinsic modes, we observed that the first IM exhibited a dominant frequency, centered around zero. Similarly, this mode showed a higher relative energy contribution among all participants. IMs 2 to 5 appeared within the neurophysiological 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 their corresponding bandwidths. We found that those results were relatively consistent among participants and experiments.

The remaining modes spanned across high frequencies. As discussed in Section 2.3, 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—appears close to the primary respiratory-related harmonic, while mode 9—with a central frequency of ~520 mHz– closely coincided with the first harmonic of the cardiac pulsations (Cordes et al., 2001; Yuen et al., 2019).

Regarding the relative energy contribution among IMs, illustrated in Figure 2b, we observed that, in general, modes with lower frequencies exhibited higher energy, while those with higher frequencies had decreasing energy contributions. Specifically, we observed that the first IM captured residual trends, characterized for having a central frequency around zero (see Table 2). This indicates that this IM is primarily capturing slow signal drifts (Huettel et al., 2009; Power et al., 2011), rather than former oscillatory components. In this way, unlike other general low-frequency components (Tong et al., 2019), the first IM primarily appeared to capture trends and slow signal drifts (Huettel et al., 2009; Power et al., 2011; Maknojia et al., 2019), which have consistently been attributed to a combination of physiological fluctuations (Huettel et al., 2009; Power et al., 2011), head motion residuals (Maknojia et al., 2019), and scanner instabilities (Smith et al., 1999). In contrast, the high-frequency components, over the neurophyisiological band, which are associated to respiratory and cardiac harmonics, showed comparatively lower energy contributions. The reason for this is that these high-frequency components are more localized (e.g., in the cerebrospinal fluid or blood vessels), resulting in a relatively smaller contribution to the overall fMRI signal, compared to other physiological components, such for example head movements, that affect the full brain.

Figure 5 depicts the average FC patterns associated only with neurophysiological modes. We obtained those FC patterns by averaging them across all participants, as they exhibited high similarities. 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 for false discovery rate adjusted to p < 0.001, as used by Romanello et al. (2022), as provide robust statistical analysis without assumptions about the underlying data distributions. For convenience, we arranged the ROIs according to the leading module (left and right), following the order reported in Supplementary Table 1.

Overall, by examining Figure 5 we observed some overlap in connectivity patterns across IM 2–5 within each experiment. To further investigate this observation, we examined the FC patterns associated with IM 1 and IM 6–10 (see Supplementary Figures 3, 4). Notably, the strong consistent connectivity patterns observed across IM 2–5 were largely absent in the higher-frequency components (IM 6–10).

The central frequency and energy distributions for MVMD in Figure 2 displayed similar trends across all experiments. However, closer examination revealed some differences between the experiments. In the gambling experiment (Figure 2), for instance, IMs 5 and 6 exhibited significantly higher central frequencies than those in other experiments.

Similary, IMs 6 and 7 also contributed significantly to the fMRI signals in the gambling experiment. These findings indicate faster physiological rhythms during the gambling experiments, which we attribute to higher arousal levels. The higher central frequency, especially in IM 6, points to faster respiratory cycles. IM 9 also showed increased energy, consistent with the first respiratory harmonic. However, when compared to the other experiments, this difference was not statistically significant.

Low-frequency IMs revealed additional differences in energy distributions across fMRI experiments. In the resting-state experiment, energy appeared concentrated primarily in the first IM, then rapidly shifted to IMs 2 and 3 before vanishing at higher frequencies. In contrast, in task-related experiments, the contribution of the first IM was considerably lower, while IMs 2-4 exhibited consistently higher energy. This pattern suggests that task-related experiments involve increased energy demands within the neurophysiological band and reflect the heightened neuronal activity required to perform cognitive tasks.

These results also align with previous research indicating that cognitive engagement leads to stronger BOLD signal changes due to localized oxygen consumption and vascular responses (Power et al., 2011). Consequently, our findings underscore the fundamental contrast between minimal cognitive engagement during resting-state scans and more complex neuronal activation patterns in task-based experiments, which ultimately leads to higher overall energy contributions within these intrinsic modes.

Regarding the FC study, Figure 5 shows the average FC patterns associated with neurophysiological IMs from MVMD. Overall, we identified some similarities between the experiments. The presence of similar connectivity patterns across multiple frequency bands likely reflects common underlying anatomical connectivity structures that support brain activity across the neurophysiological spectrum. These findings align with previous static FC studies (Bolton et al., 2020b) that the presence of similar connectivity patterns across multiple frequency bands likely reflect a common underlying connectivity structure, large-scale nets of interconnected neurons that mediates and support the signaling between brain regions (Hermundstad et al., 2013). Furthermore, the similarity of FC patterns among participants declined with increasing frequency (Supplementary Figure 4), supporting our hypothesis that higher-frequency components are more likely driven by physiological noise rather than neurophysiological processes.

On the other hand, the FC patterns associated with IM 1 presented a more complex behavior (Supplementary Figure 3), exhibiting some consistency across participants with partial overlap to IM 2–5 patterns. However, it also exhibits large individual variability as shown in Figure 6, this indicates that despite having some common patterns, it also contains large individual differences, which reflect the interplay and mixture of several components, such as scanner instabilities and motion artifacts (Maknojia et al., 2019).

Upon closer inspection, we also uncovered some differences, which provide further opportunities to explore the underlying brain mechanisms. For instance, IM 2 exhibited extensive activity across several brain networks, including the default mode network (DMN), the temporal and occipital modules. When we analyzed the resting-state patterns, we observed more significant correlations between the occipital and temporal modules than in task-related experiments, suggesting broader engagement of neural networks such as the attention network.

The motor task, on the other hand, showed a higher number of significant correlations between the occipital region and the DMN. In contrast, the gambling experiments showed fewer significant correlations. This finding suggests a specific involvement of visual control processes, aligning with FC dynamics reported in other task-related experiments (Branzi et al., 2022).

Furthermore, during the gambling experiment, we noted strong left-hemisphere correlations within the limbic module. This lateralization may reflect emotional responses to gambling tasks (Wu et al., 2022). The DMN and the occipital module also exhibited significant correlations, highlighting the interplay between intrinsic processing and visual input during this condition.

Regarding IM 3, we observed more consistent patterns across experiments. Task-related paradigms consistently showed a significant connection between the frontoparietal and temporal modules, marked by left-lateralization that weakened in the resting-state condition.

Finally, IM 4 revealed minimal differences among experiments, implying that higher-frequency neurophysiological IMs exert relatively low influence on task-related activity. We hypothesize that other paradigms designed to challenge higher cognitive functions might uncover distinct connectivity patterns at these frequencies. Future research should explore whether more demanding tasks or additional cognitive loads could alter high-frequency FC contributions.

As with any other study, there are some limitations: (1) We analyzed three fMRI datasets with similar imaging protocols. Future experiments with a broader range of datasets and conditions would help us to fully understand the applicability of the proposed approach. (2) Exploring dynamic FC aspects using MMD was beyond the scope of this study. A thorough examination of connectivity dynamics following the proposed approach and how they relate to other dynamic FC approaches could be explored in future work. (3) We focused on a single brain atlas. A comparison with other alternative brain atlases could also provide additional insights. (4) Our study focused on the ROI level. Exploring the voxel-level application of MFC could yield fine-grained insights into brain FC. (5) The current MVMD algorithms are applied individually. More advanced algorithms tailored to fMRI data and group-level strategies that integrate information from multiple participants at the same time could be investigated in future research.(6) In this study, we analyzed the three different fMRI experiments separately. As a direction of future work, we suggest performing a more comprehensive joint analysis, as this would enable more robust cross-experimental comparisons.