A Systematic Review of Data-Driven Multimodal Fusion Methods for Simultaneous EEG-fMRI

1. Introduction

Brain function emerges from coordinated neural processes that unfold across space and time, yet no single neuroimaging modality can fully capture this complexity. Electroencephalography (EEG) provides millisecond-scale temporal resolution but limited spatial specificity [1], whereas functional magnetic resonance imaging (fMRI) offers whole-brain spatial coverage with comparatively slow temporal dynamics [2]. Simultaneous EEG-fMRI acquisition provides a unique opportunity to observe these complementary aspects of brain activity within a unified framework.

EEG–fMRI fusion can be viewed as an inverse problem: recovering latent neural processes that are only partially observed through electrophysiological and hemodynamic signals. EEG reflects synchronized electrical activity of neuronal populations, while fMRI measures hemodynamic responses shaped by neurovascular coupling [3,4]. Although both signals originate from shared neural activity, they differ in spatial scale, temporal dynamics, and measurement noise. This creates a scenario in which the underlying neural processes are jointly expressed but incompletely captured in each modality [5,6]. As illustrated in Figure 1, EEG and fMRI provide complementary but incomplete observations of underlying neural processes, and fusion aims to recover their shared latent structure.

A central challenge in EEG-fMRI fusion is therefore to identify the shared latent structure that explains both electrophysiological and hemodynamic signals, while accounting for modality-specific variability and temporal misalignment. Over the past two decades, a diverse set of data-driven approaches has been developed to address this problem, ranging from linear decomposition methods such as independent component analysis (ICA) [7,8] and canonical correlation analysis (CCA) [9,10] to tensor-based models [11] and nonlinear deep learning frameworks [12]. These methods differ not only in their mathematical formulation but also in their implicit assumptions about how neural activity is represented across modalities.

Despite substantial progress, the field remains fragmented with different frameworks targeting different aspects of cross-modal coupling. This complicates comparison across studies and limits unifying interpretation. Key challenges include interpreting latent components, neurovascular effects, and spurious associations driven by noise or preprocessing. Despite these challenges, the field has grown rapidly over the past two decades. As shown in Figure 2, the number of published studies on simultaneous EEG-fMRI has increased steadily since the early 2000s, with accelerated growth in recent years.

In this review, we provide a structured and conceptually grounded overview of data-driven EEG-fMRI fusion methods developed between 2000 and 2025. We emphasize data-driven methods because of their broader applicability and reduced reliance on prior assumptions, and we focus on simultaneous acquisition because asynchronous EEG and fMRI recordings measure different physical mechanisms. Our goal is to provide a clear and structured overview of existing fusion approaches, explain how they extract and interpret shared neural information, and highlight their individual strengths and limitations.

We further examine how these methods are evaluated, the insights they have provided in cognitive and clinical neuroscience, and the limitations that continue to constrain their interpretability and generalizability. Finally, we highlight emerging directions, including nonlinear modeling, flexible coupling frameworks, and large-scale group-level fusion, that may enable EEG-fMRI integration to move beyond methodological development toward reproducible neuroscience insights and clinically meaningful biomarkers.

2. Search Strategy and Study Selection

A structured literature search was conducted to identify studies on simultaneous EEG-fMRI fusion. The search spans between approximately 2000 and 2025, a period of simultaneous EEG-fMRI acquisition and the development of modern data-driven fusion techniques. Publications were retrieved from Google Scholar and PubMed using combinations of three concept blocks: (i) EEG keywords: (“EEG" OR “electroencephalography" OR “ERP"), (ii) fMRI keywords: (“fMRI" OR “functional MRI" OR “functional magnetic resonance imaging" OR “BOLD"), and (iii) fusion keywords: (“fusion" OR “integration" OR “multimodal" OR “joint analysis" OR “coupling" OR “data-driven" OR “latent representation" OR “linked" OR “shared" OR “common space" OR “joint ICA" OR “parallel ICA" OR “mCCA" OR “tensor" OR “decomposition" OR “deep learning" OR “representation"). The initial search returned 419 records. Titles and abstracts were screened to exclude duplicate records, unimodal neuroimaging studies, non-neuroimaging multimodal studies, dissertations, preprints, and dataset only papers. Full-text review was then performed to remove trimodal integration, review or technical papers, and studies without simultaneous EEG-fMRI acquisition. In the final stage, we retained only studies employing data-driven subspace learning approaches for multimodal integration. Methods based on predefined biophysical constraints were excluded, as they fall outside the scope of this review, which focuses on learning shared latent representations directly from data. Following these criteria, $\mathbf{83}$ studies were included for detailed analysis. The full selection process is summarized in the PRISMA diagram (Figure 3).

3. Methods

Data-driven EEG–fMRI fusion identifies shared patterns across modalities with minimal prior assumptions. In this review, we grouped these methods into symmetric and asymmetric fusion [13]. Symmetric fusion methods treat EEG and fMRI equally and use joint information from both modalities. Asymmetric fusion methods typically use one modality as the prior to inform the analysis of another modality. Symmetric and asymmetric fusion further grouped into: (1) factorization, where a method aims to align latent representations of EEG and fMRI data; (2) translation, where a method aims to reconstruct EEG and fMRI data. Across both frameworks, fusion models can also be grouped as linear and nonlinear, depending on the mapping function from data to representations is linear or nonlinear. We provide a structured summary of fusion methods in Figure 4 and review each method in the following sections.

3.1. Background

3.1.1. Independent Component Analysis (ICA)

Independent component analysis (ICA) aims to separate multivariate data into a set of mutually independent sources. It assumes that the observed neural signal is a linear mixture of independent functional/neural sources that combine to form the measured data. For the observed signal, $\mathbf{X}$ can be represented as:

$$\mathbf{X}=\mathbf{A}\mathbf{S},$$

(1)

where $\mathbf{A}$ is the mixing matrix that represents the way sources are linearly combined and $\mathbf{S}$ contains the independent source signals. ICA estimates the unmixing matrix ${\mathbf{W}}^{-1}$ to recover independent sources:

$$\mathbf{S}=\mathbf{W}\mathbf{X}.$$

(2)

The fundamental assumption of ICA is that the components of $\mathbf{S}$ are typically independent with at most one Gaussian source. $\mathbf{W}$ is estimated by maximizing the independence using kurtosis, negentropy, or mutual information minimization.

EEG and fMRI signals can be viewed as mixtures of multiple neural and non-neural sources [14,15]. Applying ICA on EEG and fMRI generates separate components for each modality and then links them using correlation, regression, or source-localization approaches, has been widely used in EEG-fMRI experiments [16,17,18,19,20,21,22,23,24,25,26,27,28]. Recently, ICA has been extended to study time-resolved networks, to examine how these networks change over time and co-vary with EEG spectral power [29,30]. Group-level ICA frameworks have also been extended to EEG data, enabling estimation of components that are consistent across subjects and facilitating group inference through back-reconstruction approaches [31].

3.1.2. Canonical Correlation Analysis (CCA)

Canonical component analysis (CCA) essentially generalizes a large set of statistical tests by finding the maximally correlated coupling between two sets of examples. CCA can be thought of as finding a linear combination for each modality or subject, such that we find a shared latent representation of the data. Specifically, if we have two sets of observed signals ${\mathbf{X}}_{\mathrm{EEG}}$ and ${\mathbf{X}}_{\mathrm{fMRI}}$, CCA assumes that each can be represented as follows:

$${\mathbf{X}}_{\mathrm{EEG}}={\mathbf{A}}_{\mathrm{EEG}}\mathbf{S},{\mathbf{X}}_{\mathrm{fMRI}}={\mathbf{A}}_{\mathrm{fMRI}}\mathbf{S}.$$

(3)

Similar to ICA, $\mathbf{S}$ is the shared source between the two sets of observed signals, but ${\mathbf{A}}_{\mathrm{EEG}}$ and ${\mathbf{A}}_{\mathrm{fMRI}}$ are loading parameters that linearly map the data into the high-dimensional space of the observations. In essence, CCA assumes that data is generated from the same D-dimensional sources, and in case of EEG-fMRI fusion, these sources are assumed to be neural activity.

3.1.3. Deep Learning

Deep learning models capture nonlinear relationships between EEG and fMRI by mapping each modality into a latent space:

$${\mathbf{Z}}_{\mathrm{EEG}}=f_{{\theta }_{\mathrm{EEG}}}\left({\mathbf{X}}_{\mathrm{EEG}}\right),{\mathbf{Z}}_{\mathrm{fMRI}}=f_{{\theta }_{\mathrm{fMRI}}}\left({\mathbf{X}}_{\mathrm{fMRI}}\right),$$

(4)

where $f_{\theta }$ represents a nonlinear function parameterized by neural network weights $\theta$, and $\mathbf{Z}$ denotes the learned latent representations. These latent features play a role similar to components in ICA or CCA, but they are learned through optimization rather than predefined statistical assumptions.

3.2. Symmetric Fusion Methods

3.2.1. Joint ICA (jICA)

Joint ICA (jICA) is a symmetric, data-driven fusion framework designed to integrate EEG and fMRI [7]. jICA applies classical ICA on concatenated multimodal features [32]. In jICA, each modality is arranged as a subject-by-feature matrix as:

$${\mathbf{X}}_{\mathrm{EEG}}\in {\mathbb{R}}^{N\times M_1},{\mathbf{X}}_{\mathrm{fMRI}}\in {\mathbb{R}}^{N\times M_2},$$

(5)

where N is the number of subjects, $M_1$ is the EEG feature dimension, and $M_2$ is the fMRI feature dimension. EEG features typically consist of ERP amplitudes or spectral power measures, while fMRI features are usually voxelwise contrast maps or activation patterns. The modality matrices are variance-normalized and concatenated:

$$\mathbf{X}=\left[{\mathbf{X}}_{\mathrm{EEG}}{\mathbf{X}}_{\mathrm{fMRI}}\right]\in {\mathbb{R}}^{N\times (M_1+M_2)}.$$

(6)

jICA assumes a shared linear mixing model:

$$\mathbf{X}=\mathbf{A}\mathbf{S},$$

(7)

where $\mathbf{A}\in {\mathbb{R}}^{N\times K}$ is the shared mixing matrix and $\mathbf{S}\in {\mathbb{R}}^{K\times (M_1+M_2)}$ contains the joint sources. Partitioning the source matrix yields:

$$\mathbf{S}=\left[{\mathbf{S}}_{\mathrm{EEG}}{\mathbf{S}}_{\mathrm{fMRI}}\right],$$

(8)

linking EEG and fMRI patterns for each independent component.

jICA has been extensively applied in both cognitive and clinical paradigms. It has been widely used in estimating joint ERP-BOLD components in visual cognitive tasks [33,34], auditory oddball experiments [7,35,36], realistic decision-making paradigms [37], cognitive control processes [38], and clinical populations such as schizophrenia [39]. Moreover, simulation studies have validated the effectiveness of jICA for multimodal joint decomposition under controlled conditions [32]. Overall, jICA produces interpretable paired spatiotemporal components with shared subject loadings, making it particularly suitable for identifying cross-modal patterns in task-based EEG-fMRI data. However, the method assumes a single shared mixing matrix across modalities and relies on linear mixing assumptions. These constraints, together with sensitivity to feature scaling, can limit its flexibility in modeling partially shared or modality-specific processes.

3.2.2. Group Joint ICA (gjICA)

To move beyond single subject analyses, group joint ICA (gjICA) was developed as a group-level framework that estimates multimodal components at the group level followed by estimation of subject-level components that correspond across subjects [40]. In gjICA, subject-level EEG and fMRI data are first combined and reduced using principal component analysis (PCA) and then concatenated across subjects before performing group ICA. Let ${\mathbf{X}}_{\mathrm{EEG}}^{\left(i\right)}$ and ${\mathbf{X}}_{\mathrm{fMRI}}^{\left(i\right)}$ denote the feature matrices of subject i. For each subject, the modalities are first combined into a joint representation:

$${\mathbf{X}}^i=\left[{\mathbf{X}}_{\mathrm{EEG}}^{\left(i\right)}{\mathbf{X}}_{\mathrm{fMRI}}^{\left(i\right)}\right].$$

(9)

After dimensionality reduction using PCA, the reduced subject data ${\tilde{\mathbf{X}}}_i$ are concatenated across all N subjects to form the group data matrix:

$${\mathbf{X}}_G=\left[\begin{array}{c}{\tilde{\mathbf{X}}}^{\left(1\right)}\\ {\tilde{\mathbf{X}}}^{\left(2\right)}\\ ⋮\\ {\tilde{\mathbf{X}}}^{\left(N\right)}\end{array}\right].$$

(10)

jICA then decomposes the concatenated matrix as

$${\mathbf{X}}_G={\mathbf{A}}_G{\mathbf{S}}_G,$$

(11)

where ${\mathbf{S}}_G$ represents the joint independent components shared across subjects and modalities, and ${\mathbf{A}}_G$ contains the subject-specific loading parameters describing the contribution of each component to each subject.

This decomposition produces a common set of joint components that can be interpreted at the group level, making it possible to identify shared multimodal brain patterns across the entire sample and to perform statistical comparisons between groups. At the same time, gjICA enables estimation of subject-specific expressions of the joint components through the loading matrix ${\mathbf{A}}_G$. Furthermore, gjICA introduces joint time courses that are shared across EEG and fMRI representations, providing a unified temporal profile associated with the coupled multimodal patterns.

3.2.3. Parallel ICA

Parallel ICA [41,42] extends standard ICA by jointly optimizing statistical independence within each modality and correlation across modalities. Unlike jICA, parallel ICA does not assume a shared mixing matrix. Instead, EEG and fMRI are decomposed separately, while a correlation term is introduced to maximize the association between selected component loadings across modalities. Formally, let ${\mathbf{A}}_{\mathsf{fMRI}}$ and ${\mathbf{A}}_{\mathsf{EEG}}$ be mixing matrices within fMRI and EEG, respectively. The correlation between a single column i in the fMRI and a single column j in the EEG is defined as

$$\mathsf{Corr}({\mathbf{A}}_{\mathsf{fMRI},i},{\mathbf{A}}_{\mathsf{EEG},j})={\displaystyle \frac{\mathsf{Cov}({\mathbf{A}}_{\mathsf{fMRI},i},{\mathbf{A}}_{\mathsf{EEG},j})}{\mathsf{Std}\left({\mathbf{A}}_{\mathsf{fMRI},i}\right)·\mathsf{Std}\left({\mathbf{A}}_{\mathsf{EEG},j}\right)}}.$$

(12)

The optimization objective combines entropy maximization for ICA in each modality with the cross-modal correlation constraint:

$$\mathcal{L}=H\left({\mathbf{Y}}_{\mathsf{EEG}}\right)+H\left({\mathbf{Y}}_{\mathsf{fMRI}}\right)+\lambda \mathsf{Corr}({\mathbf{A}}_{\mathsf{fMRI},i},{\mathbf{A}}_{\mathsf{EEG},j}),$$

(13)

where $\lambda$ controls the tradeoff between within-modality independence and cross-modality coupling. The ICA outputs are typically written as

$${\mathbf{Y}}_m=g\left({\mathbf{U}}_m\right),{\mathbf{U}}_m={\mathbf{W}}_m{\mathbf{X}}_m+{\mathbf{W}}_{0,m},m\in \{\mathsf{EEG},\mathsf{fMRI}\},$$

(14)

with ${\mathbf{W}}_m$ the unmixing matrix, ${\mathbf{W}}_{0,m}$ a bias term, and $g(·)$ a nonlinear activation function, often logistic:

$$g\left(u\right)={\displaystyle \frac{1}{1+e^{-u}}}.$$

(15)

In practice, the entropy terms in Eq. 13 are optimized using natural gradient learning [43], while the correlation term is optimized separately, for example, using steepest descent with step sizes determined by Wolfe conditions [44]. To prevent overfitting to a single objective term, the learning rates are typically annealed during training.

Parallel ICA was introduced for multimodal neuroimaging fusion by Liu et al. [41] and further formalized in Pearlson et al. [42]. In cognitive neuroscience, parallel ICA has been used in auditory target detection tasks [8]. Simulation-based validation studies demonstrated its robustness under controlled multimodal conditions [45]. Beyond task paradigms, parallel ICA has been applied to resting-state data to investigate neurophysiological correlates of intrinsic networks [46]. In clinical applications, it has been used in epilepsy to explore coupling between electrophysiological signatures and hemodynamic networks [47]. Parallel ICA provides greater flexibility than jICA because it does not assume a single shared mixing matrix. Each modality retains its own independent components, while cross-modal relationships are enforced at the level of subject loadings. This makes it well suited for partially shared neural processes.

3.2.4. Bidirectional Independent Component Averaged Representation (BICAR)

Unlike parallel ICA, which links components through correlation constraints during the decomposition process, bidirectional independent component averaged representation (BICAR) identifies cross-modal correspondences after independent decompositions and stabilizes the results through averaging across multiple ICA runs [48]. In BICAR, ICA is applied repeatedly to each modality independently. For each run k, independent decompositions are performed:

$${\mathbf{X}}_{\mathsf{EEG}}={\mathbf{A}}_{\mathsf{EEG}}^{\left(k\right)}{\mathbf{S}}_{\mathsf{EEG}}^{\left(k\right)},{\mathbf{X}}_{\mathsf{fMRI}}={\mathbf{A}}_{\mathsf{fMRI}}^{\left(k\right)}{\mathbf{S}}_{\mathsf{fMRI}}^{\left(k\right)}.$$

(16)

Cross-modal correspondences are then identified by comparing EEG components with fMRI mixing profiles. To improve robustness, matched component pairs across multiple ICA runs are averaged to obtain stable multimodal components:

$$\overline{\mathbf{S}}={\displaystyle \frac{1}{K}}{\displaystyle \sum _{k=1}^K}{\mathbf{S}}^{\left(k\right)},$$

(17)

where $\overline{\mathbf{S}}$ represents the averaged component estimate. This averaging step reduces variability due to random initialization and improves reproducibility of the multimodal decomposition. In simultaneous EEG–fMRI fusion, BICAR has successfully identified biologically meaningful components related to visual processing, motor activity, and attentional networks [49].

3.2.5. Connectivity ICA (connICA)

Connectivity ICA (connICA) is a data-driven framework designed to extract hybrid connectivity patterns jointly expressed across multiple modalities [50]. Instead of decomposing raw signals or activation maps, connICA operates on functional connectivity matrices derived from each modality. In the context of EEG-fMRI fusion, whole-brain functional connectivity (FC) matrices are first constructed for both modalities [51]. For a subject i, the connectivity matrices can be written as:

$${\mathbf{C}}_{\mathsf{fMRI}}^{\left(i\right)}=\mathrm{corr}\left({\mathbf{T}}_{\mathsf{fMRI}}^{\left(i\right)}\right),{\mathbf{C}}_{\mathsf{EEG}}^{\left(i\right)}=\mathrm{conn}\left({\mathbf{T}}_{\mathsf{EEG}}^{\left(i\right)}\right),$$

(18)

where $\mathbf{T}$ denotes the regional time series and $\mathrm{conn}(·)$ represents an EEG connectivity measure. Each connectivity matrix is then vectorized by extracting the upper triangular elements:

$${\mathbf{x}}^{\left(i\right)}=\left[\mathrm{vec}\left({\mathbf{C}}_{\mathsf{fMRI}}^{\left(i\right)}\right),\mathrm{vec}\left({\mathbf{C}}_{\mathsf{EEG}}^{\left(i\right)}\right)\right].$$

(19)

The subject-level vectors are stacked across all N subjects to form a data matrix:

$$\mathbf{X}=\left[\begin{array}{c}{\mathbf{x}}^{\left(1\right)}\\ {\mathbf{x}}^{\left(2\right)}\\ ⋮\\ {\mathbf{x}}^{\left(N\right)}\end{array}\right].$$

(20)

Dimensionality reduction is first performed using PCA and ICA is then applied to the reduced representation:

$${\mathbf{U}}^{\top }\mathbf{X}=\mathbf{A}\mathbf{S},$$

(21)

where $\mathbf{U}$ contains the principal components and $\mathbf{S}$ contains independent connectivity components (connectivity traits) and $\mathbf{A}$ represents subject-specific mixing weights that quantify how strongly each subject expresses each hybrid connectivity pattern. Each component thus corresponds to a paired EEG-fMRI connectivity structure together with subject loadings that capture inter-subject variability in multimodal network organization.

3.2.6. Multi-Set CCA (mCCA)

Multi-set CCA (mCCA) extends canonical correlation analysis from two datasets to multiple datasets, making it well-suited for group-level EEG-fMRI fusion [9,10]. In simultaneous EEG-fMRI studies with N subjects, each participant contributes two datasets (EEG and fMRI), resulting in $2N$ datasets that can be analyzed jointly within a single framework. This extension allows identification of neural processes that are shared across modalities and consistent across subjects. Let ${\mathbf{X}}_i\in {\mathbb{R}}^{T\times p_i}$ denote dataset i for $i=1,\dots ,2N$, where T is the number of samples and $p_i$ is the number of features. mCCA seeks projection vectors ${\mathbf{w}}_i$ such that the canonical variates:

$${\mathbf{y}}_i={\mathbf{X}}_i{\mathbf{w}}_i$$

(22)

are maximally correlated across all datasets. The optimization problem can be written as:

$$\underset{{\mathbf{w}}_1,\dots ,{\mathbf{w}}_{2N}}{max}\sum _{i\ne j}\mathrm{corr}({\mathbf{X}}_i{\mathbf{w}}_i,{\mathbf{X}}_j{\mathbf{w}}_j),$$

(23)

subject to normalization constraints on the projected variables.

In EEG-fMRI fusion, the resulting canonical variates represent shared patterns of variation that are consistent across subjects and modalities. This enables identification of neural processes that jointly explain variability in EEG and fMRI while naturally accounting for subject-level differences.

3.2.7. Tensor Decomposition-Based Fusion

Tensor decomposition provides a principled framework for multimodal fusion by preserving the multi-way organization of neuroimaging data [52]. EEG naturally forms higher-order tensors (e.g., subjects × channels × time or frequency), and fMRI can be organized as a subject × space matrix. Central to the tensor decomposition family is the canonical polyadic decomposition (CPD) or PARAFAC [53,54,55,56,57], which expresses an N-way (multi-dimensional) tensor as a sum of rank-one outer products across dimensions:

$$\mathcal{X}\approx {\displaystyle \sum _{r=1}^R}{\lambda }_r{\mathbf{a}}_r^{\left(1\right)}\circ {\mathbf{a}}_r^{\left(2\right)}\circ \dots \circ {\mathbf{a}}_r^{\left(N\right)}.$$

(24)

Because PARAFAC assumes approximate multilinearity, separability across dimensions, and additive noise [58], it is well suited for EEG–fMRI fusion, where shared neural generators may manifest differently across modalities but can still be represented through common latent factors. A widely-used fusion strategy in this context is coupled matrix–tensor factorization (CMTF), in which EEG data are modeled as a tensor and fMRI data as a matrix with shared factor matrices across one or more modes [59]. Let the EEG data be represented as:

$${\mathcal{X}}_{\mathrm{EEG}}\in {\mathbb{R}}^{I\times J\times K},$$

(25)

with the PARAFAC model:

$${\mathcal{X}}_{\mathrm{EEG}}\approx {\displaystyle \sum _{r=1}^R}{\mathbf{s}}_r\circ {\mathbf{c}}_r\circ {\mathbf{t}}_r=〚\mathbf{S},\mathbf{C},\mathbf{T}〛.$$

(26)

Similarly, fMRI data are arranged as a subject × voxel matrix X,

$${\mathbf{X}}_{\mathrm{fMRI}}\in {\mathbb{R}}^{I\times V},$$

(27)

are modeled via:

$${\mathbf{X}}_{\mathrm{fMRI}}\approx \mathbf{S}{\mathbf{A}}^{\top },$$

(28)

where the subject-mode factor matrix $\mathbf{S}$ is shared across EEG and fMRI, ensuring that each component jointly captures: (1) EEG spatial scalp patterns ($C_r$), (2) EEG temporal signatures ($t_r$), (3) fMRI spatial maps ($a_r$), and (4) subject loadings ($S_r$). Joint estimation typically minimizes:

$$\underset{\mathbf{S},\mathbf{C},\mathbf{T},\mathbf{A}}{min}{∥{\mathcal{X}}_{\mathrm{EEG}}-〚\mathbf{S},\mathbf{C},\mathbf{T}〛∥}_F^2+\alpha {∥{\mathbf{X}}_{\mathrm{fMRI}}-\mathbf{S}{\mathbf{A}}^{\top }∥}_F^2,$$

(29)

with $\alpha$ controlling the balance between modalities. Extensions include PARAFAC2 for handling time variability in EEG [60], soft-coupled tensor decompositions [61,62,63] allowing flexible alignment [62,63], and correlation– [64,65] or mutual-information–guided coupling for more adaptive matching between modalities [66,67]. PARAFAC-based tensor fusion has been applied extensively to simultaneous EEG–fMRI experiments to extract unified multimodal components.

3.2.8. Deep Learning-Based Symmetric Fusion

Symmetric fusion strategies based on deep learning commonly learn shared representations of EEG and fMRI through multimodal encoders, latent alignment, attention mechanisms, or graph-based integration. In these models, both modalities contribute to a common representation space, enabling the network to capture complementary temporal information from EEG and spatial information from fMRI. In EEG-fMRI fusion, the objective function typically aims to maximize similarity between latent features ${\mathbf{Z}}_{\mathrm{EEG}}$ and ${\mathbf{Z}}_{\mathrm{fMRI}}$ extracted from the two modalities:

$$max\mathcal{L}=max\mathrm{sim}\left({\mathbf{Z}}_{\mathrm{EEG}},{\mathbf{Z}}_{\mathrm{fMRI}}\right)=\underset{\theta ,\varphi }{max}\mathrm{sim}\left({\mathbf{f}}_{\theta }\left({\mathbf{X}}_{\mathrm{EEG}}\right),{\mathbf{g}}_{\varphi }\left({\mathbf{X}}_{\mathrm{fMRI}}\right)\right).$$

(30)

In practice, this similarity objective is often combined with task-specific classification or reconstruction losses:

$${\mathcal{L}}_{\mathrm{total}}={\mathcal{L}}_{\mathrm{task}}+\lambda {\mathcal{L}}_{\mathrm{align}},$$

(31)

where ${\mathcal{L}}_{\mathrm{align}}$ encourages agreement between the learned EEG and fMRI embeddings.

Multimodal decoding enhances visual object categorization by leveraging complementary temporal and spatial information [12], while graph-based fusion of EEG and fMRI connectivity improves prediction of antidepressant response in major depressive disorder [68]. In inner speech recognition, cross-perception architectures combining modality-specific encoders, shared representations, and attention-based fusion outperform unimodal baselines [69]. Multimodal embedding frameworks further extend these approaches by learning unified latent representations across neural modalities. For example, BrainFLORA captures consistent concept-level structure across EEG, MEG, and fMRI [70]. Deep learning has also enabled symmetric translation-based fusion, where hierarchical transcoding models cyclically map EEG and fMRI to recover high spatiotemporal latent representations [71,72]. Overall, deep learning-based symmetric fusion provides a flexible framework for modeling nonlinear and partially shared EEG–fMRI relationships. Compared with linear approaches such as ICA and CCA, these methods capture more complex cross-modal structure but require larger datasets and raise challenges in interpretability, generalization, and biological validity.

3.3. Asymmetric Fusion

Asymmetric fusion uses one modality to constrain or inform the other. Unlike symmetric fusion, asymmetric fusion often relies on explicit (biophysical) assumptions about how EEG and fMRI are related. In contrast, symmetric fusion approaches often operate under relatively lenient assumptions. In general, there are three main fields within asymmetric fusion. First, the use of identifiable temporal events or EEG traces to perform asymmetric fusion. Second, there are methods that define specific resting-state features in one signal and relate them to the other signal. Lastly, there are translation-based methods, where one modality is used to reconstruct the raw signal from the other, without any pre-defined features.

3.3.1. Event-Based Asymmetric Fusion

Many asymmetric fusion studies use either derived events or task-based designs to fuse fMRI with EEG data. One of the two modalities is then often regressed to these events, and asymmetrically integrated with the other modality. There are two broad types of events that are used for asymmetric fusion, namely ictal events, and task-based events. The timing of the ictal events is based on a clinical assessment of the seizure onset, whereas the task-based events are often specific block designs.

Ictal Events: Ictal events, often referred to as interictal epileptic discharges (IEDs), are short abnormal EEG waveforms that signal seizure onset. Because these events are temporally distinct, they are well suited for EEG–fMRI fusion: IEDs are detected in EEG, and their timings are used to examine corresponding BOLD responses in fMRI. Lee et al. [73] applied ICA to detect IED-related EEG components and used these signals to define hemodynamic response functions (HRFs) in a GLM to localize seizure-related fMRI activity. LeVan et al. [74] instead applied spatial ICA to fMRI and modeled EEG-derived ictal events as boxcar regressors, estimating HRF timing separately for each region. To further relax HRF assumptions, Masterton et al. [75] introduced event-related ICA (eICA), applying spatial ICA to deconvolved fMRI signals linked to IED timings and demonstrating improved identification of seizure onset regions. To automate event detection, Omidvarnia et al. [76] proposed a matched filtering method that simplified EEG–fMRI fusion and, in some cases, outperformed manual IED identification.

Task-Based Events: Task-based fusion similarly extracts EEG-derived neural signatures that are regressed to fMRI to localize task-related activity. Ahmad et al. [77] averaged EEG responses to visual stimuli and used GLM regression with a canonical HRF to relate EEG features to fMRI activations, showing improved classification when combining modalities. Muller et al. [78] used ERP components such as P1 as regressors and localized activity to the amygdala and lateral posterior cortex. Pisauro et al. [79] constructed EEG-informed regressors reflecting trial-by-trial evidence accumulation and identified the posterior medial frontal cortex as a key region. He et al. [80] related alpha-band power to fMRI responses during gesture–speech integration using window-based analyses to track spatiotemporal coupling. Several studies extend this idea through multiple asymmetric fusion steps. Dong et al. [81,82] applied spatial ICA to fMRI and temporal ICA to EEG, matched components using the maximal information coefficient, and constrained the EEG inverse problem using task activation maps to estimate spatiotemporal components. A related approach [83] proposed a variational Bayesian framework for solving the fMRI-constrained EEG inverse problem. In this method, EEG and fMRI are asymmetrically fused multiple times during inference and can be applied to both task-based and resting-state data.

3.3.2. Resting-State Asymmetric Fusion

Event-based asymmetric fusion is powerful but limited to datasets with clearly defined events. Many resting-state applications, including psychiatric and neurological biomarker development, lack such event markers. To address this, Ferdowsi et al. [60] proposed PARAFAC2 for resting state data. They used Rolandic beta band power fluctuations, convolved with a canonical HRF, and demonstrated that their method can accurately identify brain regions associated with post-movement rebound, without requiring task timings. Similarly, Feige et al. [84] regressed multiple EEG power bands against deconvolved fMRI signals at varying time lags and found that peaks in deconvolved fMRI activity often precede increases in EEG power across several frequency bands, revealing spatial patterns consistent with functional segregation in the brain. Double asymmetric fusion approaches have also been extended to resting-state analysis. For example, Lei et al. [85] proposed a parallel empirical Bayes framework that applies temporal ICA to EEG and spatial ICA to fMRI, then links the resulting components in both spatial and temporal domains. This approach solves an fMRI-constrained EEG inverse problem while estimating HRF responses from EEG temporal features. Although validated using simulation data, the framework was designed with resting-state applications in mind. Finally, [86] combined source-reconstructed EEG connectivity measures with fMRI features to predict structural connectivity derived from diffusion MRI. They showed that incorporating EEG metrics such as imaginary coherence improved prediction beyond fMRI-only models, highlighting the complementary value of multimodal fusion for studying structure–function relationships.

3.3.3. Translation-Based Asymmetric Fusion

Cross-modality translation aims to reconstruct one brain imaging modality (e.g., fMRI) from another (e.g., EEG) by learning a nonlinear mapping that preserves shared neural information [87]. This framework enables nonlinear alignment between EEG and fMRI without relying on explicit analytical models such as CCA. Two main variations of this framework exist:

Deterministic Encoder–Decoder Translation Models: Many EEG-to-fMRI (or fMRI-to-EEG) methods use deterministic encoder-decoder architectures to directly map between modalities. Convolutional neural networks (CNNs) have been widely used to learn spatial and temporal mappings and relationships between modalities. For example, Liu et al. [88] presents CNN-based transcoders, trained on simulated data, showing that deep models can learn nonlinear cross-modal relationships between EEG and fMRI. Further, Sasikumar et al. [89] designed a multi-head autoencoder model for decoding subcortical brain activity from EEG-fMRI data. Transformer-based approaches such as NeuroBolt [90] use self-attention mechanisms to capture long-range temporal dependencies in EEG and synthesize high-fidelity fMRI volumes, outperforming many baseline methods.

Probabilistic and Generative Translation Models: Another class of methods uses probabilistic and generative models that learn the distribution of one modality conditioned on the other. These include variational autoencoders (VAEs), generative adversarial networks (GANs), and diffusion probabilistic models. For example, Calhas and Henriques [91] represent one of the earliest attempts to apply VAE and $\beta$-VAE to the EEG-to-fMRI translation problem. Cheng et al. [92] proposes a GAN-based architecture for synthesizing fMRI from EEG using an adversarial training strategy that encourages the generator to produce realistic and plausible fMRI data. The model incorporates additional losses, including the cycle-consistency and a non-adversarial perception loss, to improve the quality of the generated data. Yao et al. [93] introduces a diffusion framework that conditions the reverse-diffusion process on EEG signals to synthesize fMRI and learn a shared representation.

A summary of the main characteristics, assumptions, and limitations of the reviewed fusion methods is provided in Table 1.

3.4. Other Fusion Methods

In addition to the major factorization and translation frameworks described above, several studies have proposed hybrid or specialized fusion strategies that combine multiple analytical techniques. For example, consecutive independence and correlation transform (CICT) integrates ICA and CCA to identify multimodal components while relaxing some of the constraints of joint ICA approaches [97]. Similarly, joint common spatial pattern methods have been proposed to extract discriminative multimodal sources across experimental conditions by combining EEG spatial filtering with fMRI features [98]. Recent work has also explored joint source decomposition techniques to identify shared EEG–fMRI components across subject groups [99].

Yu et al. [1] constructed multimodal brain graphs by integrating EEG spectral features with fMRI connectivity to investigate brain network organization. Other studies examine relationships between electrophysiological activity and hemodynamic signals using spectral and spatiospectral heuristic models [100,101].

Additional hybrid approaches focus on specific neurophysiological phenomena. For example, simultaneous EEG–fMRI has been used to investigate electrophysiological signatures of the fMRI global signal [102], the relationship between alpha synchrony and BOLD responses [103], and the coupling between EEG activity and dynamic fMRI connectivity using wavelet-based analyses [104]. In clinical applications, hybrid frameworks combining EEG source imaging with EEG-correlated fMRI have been used to improve localization of epileptic networks and spike-related activity [105,106]. Together, these studies illustrate a diverse set of hybrid EEG–fMRI fusion strategies that integrate decomposition, connectivity analysis, source imaging, and biophysically motivated modeling to capture complementary aspects of neural and hemodynamic brain activity.

4. Evaluation and Interpretation of EEG-fMRI Fusion

Interpreting EEG-fMRI fusion results is often more challenging than running the models themselves. In unimodal studies, EEG is typically visualized as waveforms, spectra, or topographies, and fMRI as activation maps or networks. In contrast, fusion methods produce joint representations across different spatial and temporal scales. As a result, there is no single universal evaluation metric applicable to all fusion methods. Below we summarize commonly-used evaluation strategies.

Cross-Modal Coupling and Consistency: A common approach is quantifying statistical coupling between modalities. This includes correlations between EEG-derived features, such as band-limited power, event-related potentials, or component time courses, and fMRI signals at the voxel or network level [8,16,17,18]. In symmetric fusion frameworks, associations are frequently assessed through correlations between subject-level component loadings [7,10,41]. Temporal coupling is often examined by relating EEG time courses to fMRI-derived signals while accounting for hemodynamic delays [16]. Although these measures capture statistical dependence, they do not necessarily imply shared neural generators, as similar relationships may arise from confounds, preprocessing choices, or indirect coupling mediated by neurovascular dynamics.

Spatial and Network-Level Validation: For fMRI-derived components, spatial maps are typically compared with known functional networks using measures of spatial overlap, template matching, or reproducibility across subjects and datasets [34,39,40,45]. EEG components are evaluated based on their spectral characteristics, including dominance in canonical frequency bands, as well as their scalp topographies or source-reconstructed spatial distributions [10,40,46]. In multimodal settings, a key validation step is assessing whether EEG and fMRI features jointly reflect coherent and interpretable brain networks, often by examining consistency between spectral signatures and spatial organization.

Connectivity-Based Evaluation: Connectivity analyses provide interactions between brain regions. In fMRI, functional network connectivity (FNC) is commonly-used to quantify correlations between network time courses [49,50], while EEG connectivity is assessed using measures such as coherence, phase-locking value (PLV), imaginary coherence, and weighted phase lag index (wPLI), which reduce the influence of volume conduction [51,66]. Recent work has shown that dynamic functional network connectivity (dFNC) states derived from fMRI are associated with distinct EEG spectral signatures, providing evidence that time-varying connectivity patterns reflect underlying electrophysiological processes [107,108]. Multimodal extensions, including joint functional network connectivity (jFNC), integrate these measures to capture cross-modal network interactions [40]. Statistical evaluation of connectivity patterns often involves group-level comparisons using t-tests, analysis of variance, and multiple comparison correction methods such as false discovery rate (FDR) [40,107,108]. These approaches help determine whether fusion-derived connectivity patterns reflect meaningful large-scale network organization and whether they differentiate between conditions or populations.

Reconstruction and Predictive Performance: For translation-based and deep learning models, evaluation frequently focuses on reconstruction accuracy and predictive performance. Common metrics include mean squared error (MSE), root mean squared error (RMSE), structural similarity index (SSIM), peak signal-to-noise ratio (PSNR), and correlation coefficients between predicted and observed signals [12]. In classification or regression tasks, performance is assessed using metrics such as accuracy, sensitivity, specificity, area under the receiver operating characteristic curve (AUC), and mean absolute error (MAE) [74,75]. While these metrics are useful for quantifying model performance, they primarily reflect predictive capability and do not necessarily indicate that the learned representations correspond to biologically meaningful neural processes.

Model Stability and Confound Control: Reliable interpretation also requires verifying that the results are robust and not driven by noise or confounding factors. For example, ICA-based fusion methods often assess component stability using resampling-based approaches such as ICASSO, which evaluates the reproducibility of estimated components across multiple ICA runs [109]. Statistical significance of fusion results can be evaluated using permutation testing or bootstrap resampling [110]. In addition, potential confounds such as head motion, physiological noise, or scanner artifacts should be carefully controlled during preprocessing or regression [111]. Finally, cross-validation and testing on independent datasets are commonly-used to evaluate the generalization of predictive or deep learning fusion models [112].

External and Clinical Validation: External validation assesses whether fusion-derived features relate to behavioral, cognitive, or clinical variables. This includes correlating subject-level component loadings or connectivity measures with task performance, symptom severity, or treatment outcomes, as well as evaluating group differences using statistical models that incorporate relevant covariates [7,13,36]. Such analyses are particularly important in clinical applications, where the goal is to identify robust and generalizable biomarkers.

Back-Reconstruction of Subject and Modality-Specific Components: Many fusion methods estimate components at the group level and subsequently reconstruct subject-specific representations through a process known as back-reconstruction. This step is essential for interpreting multimodal components at the individual subject level and for performing statistical comparisons across groups. Several back-reconstruction strategies have been proposed, including direct group ICA approaches (GICA1–GICA3) and indirect approaches such as spatiotemporal regression (STR) [113,114].

5. Simultaneous EEG-fMRI Datasets

The development of multimodal EEG-fMRI fusion methods has been supported by the growing availability of simultaneous EEG-fMRI datasets. Beyond a few early releases, several simultaneous EEG–fMRI datasets are now available through platforms such as OpenNeuro1. These datasets span a variety of experimental paradigms, including resting-state recordings, naturalistic audiovisual stimulation, sleep studies, cognitive and perceptual tasks, and clinical or patient populations. Most datasets include raw EEG, structural and functional MRI, and often additional physiological or behavioral measurements, typically organized in standardized formats such as BIDS.

6. Applications

In this section, we summarize key insights gained from EEG-fMRI fusion across major application areas, including neurological disorders, psychiatric disorders, and cognitive task studies.

6.1. Neurological Disorders

EEG-fMRI fusion has shown that epileptic activity is not confined to focal regions but involves distributed large-scale networks. Across studies, BOLD responses associated with interictal discharges extend beyond the presumed epileptogenic zone into cortical and subcortical systems [74,106]. Brain components with early hemodynamic responses are spatially aligned with seizure onset regions, whereas later components reflect propagation through distributed networks [74], enabling separation of onset and spread dynamics. In generalized epilepsy, multimodal analyses identified coordinated involvement of the default mode, basal ganglia, salience, and frontal networks, with both linear and nonlinear relationships to epileptic discharges [82]. These findings suggest that subcortical–cortical interactions play a central role in modulating seizure activity. Finally, integrating EEG spectral and time–frequency features with fMRI has shown that oscillatory activity maps onto distributed epileptic networks rather than isolated regions, including thalamocortical and default mode circuits, highlighting the network-level organization of epileptic dynamics [47,104].

6.2. Psychiatric Disorders

A significant amount of recent work using simultaneous EEG-fMRI in major depressive disorder (MDD) has focused on multimodal neurofeedback [115,116,117], as well as on treatment response [118] and remission [119]. These studies demonstrate that frontal EEG asymmetry and oscillatory dynamics are coupled with BOLD activity in the amygdala, reflecting impaired top–down regulation of emotional processing in MDD [115]. EEG-informed fMRI analyses further show that band-limited spectral power, ERP component (P300) can predict treatment response and remission by modulating large-scale functional networks [118,119]. Recent multimodal fusion work [68] identifies treatment-related network structure across the frontoparietal control, dorsal and ventral attention, and limbic systems, with overlapping contributions from the posterior cingulate cortex and superior temporal gyrus. The multimodal model predicts HAMD17 change more accurately than EEG- or fMRI-only models, indicating that cross-modal connectivity captures variance in treatment response not present in either modality alone.

In contrast to MDD, schizophrenia studies identified reduced ERP amplitudes and altered oscillatory activity are linked with BOLD changes in auditory cortex, frontal regions, and default mode networks, forming joint components that differentiate patients from controls and vary with symptom severity [39]. Multimodal decomposition further separates components reflecting shared EEG–fMRI structure from those dominated by a single modality, indicating the coexistence of cross-modal disruptions and modality-specific alterations [97].

6.3. Cognitive Tasks

EEG-fMRI fusion methods have been applied to identify brain signals and regions associated with cognitive functions, including auditory and visual processing, decision making, and memory encoding.

6.3.1. Auditory Processing

In the auditory oddball paradigm, Eichele et al. [8] revealed a previously undetected component reflecting auditory onset and low-level change detection, linking N1-related EEG activity with bilateral temporal and anterior cingulate fMRI response. mCCA revealed linear decreasing trends in amplitude modulations for both EEG and fMRI, with spatial activations primarily in motor, frontal, temporal, inferior parietal, and orbito-frontal areas [10]. Tensor-based fusion revealed shared components between EEG and fMRI in the auditory oddball task [66]. More recently, a mutual information-based tensor decomposition identified three shared components in the auditory oddball paradigm, exhibiting dominant frequency responses in the alpha and theta bands and linked to widespread fMRI activations across attention- and cognition-related brain regions [67].

6.3.2. Visual Processing

Similarly, the visual oddball paradigm evaluates the detection of infrequent targets among standard stimuli. Simultaneous EEG–fMRI showed a temporal cascade, with early right-lateralized frontal and occipital activations linked to attentional orienting, followed by post-response engagement of the precuneus and posterior cingulate cortex [120]. Tensor-based fusion further identified frontal, parietal, cingulate, and dorsolateral prefrontal activation linked with increased delta and theta EEG activity [62]. In a spatial attention task, contralateral visual cortex activity covaried positively with occipital gamma power, ipsilateral regions covaried inversely with alpha oscillations, and the pulvinar nucleus of the thalamus covaried with both patterns, suggesting its regulatory role in modulating visual cortex excitability [121]. These findings highlight coordinated cortical–subcortical mechanisms supporting spatial attention.

6.3.3. Decision Making

Using simultaneous EEG-fMRI, EEG signatures of evidence accumulation were linked to activity in the posterior–medial frontal cortex, which tracked trial-by-trial variability and showed coupling with the ventromedial prefrontal cortex and striatum [79]. In a complex intertemporal task, jICA identified components reflecting reward, cost, and outcome uncertainty, with corresponding EEG-fMRI patterns capturing their spatiotemporal dynamics [37].

6.3.4. Memory Encoding

Using simultaneous EEG–fMRI, long-term memory encoding is characterized by decreased beta and increased theta power in EEG, as well as strong activation in the left inferior prefrontal cortex and the right parahippocampal gyrus in fMRI, suggesting distinct roles of beta oscillations in the left inferior prefrontal cortex in semantic encoding and theta oscillations in the medial temporal lobe in item–context binding [27]. Complementing these findings, fused EEG–fMRI components were dominated by frontal theta and posterior alpha activity, with corresponding fMRI components primarily reflecting the attention network and, to a lesser extent, the default mode network in an N-back working memory task [122].

7. Discussion

7.1. From Multimodal Integration to Latent Neural Inference

A key advance in EEG–fMRI research has been the shift from modality-specific analysis to modeling shared neural processes. Early work relied on EEG-informed fMRI or post hoc correlations, whereas data-driven fusion reframed the problem as identifying latent neural activity expressed across modalities. [10,16,23,87,123]. ICA-based models enabled recovery of meaningful networks without strong priors, and extensions such as jICA and parallel ICA linked EEG and fMRI within a common framework. [7,8,14,85,123,124]. Later, CCA and mCCA reframed the problem, making it easier to study relationships across subjects [9,10]. Tensor-based methods went a step further by preserving the multi-dimensional structure of EEG data, allowing space, time, and frequency information to be modeled jointly while linking to fMRI spatial patterns [11,63,66,94]. More recently, deep learning approaches have relaxed many of the constraints of earlier models. By allowing nonlinear mappings and cross-modal translation, these methods can capture more complex relationships between EEG and fMRI and improve performance in predictive settings [12,68,70,71]. At the same time, these shift the focus toward learning representations directly from data rather than relying on predefined assumptions. Overall, these developments make EEG-fMRI fusion particularly powerful for studying complex brain function.

7.2. What Do Fusion Methods Actually Recover?

Symmetric fusion treats EEG and fMRI as equally informative and uncovers shared information present in both modalities. These approaches are useful because they do not require choosing one modality as the “driver” and can reveal multimodal patterns that would be missed by unimodal pipelines [10,123]. Within symmetric methods, jICA, parallel ICA, and mCCA-based methods explicitly maximize shared latent structure using statistical constraints such as independence and correlation [7,10,85]. Tensor-based approaches extend this idea by preserving the multi-dimensional structure of EEG data while linking it to fMRI representations, often improving interpretability when working with time–frequency features [59,62,63].

An intuitive example of what these methods recover is provided by jICA (Figure 5). Across subjects, variability in ERP waveforms and fMRI activation patterns reflects a shared neural process, which jICA decomposes into paired components linking temporal EEG dynamics with spatial fMRI maps. Recombining these components reveals a temporal sequence of spatially distinct brain activations, for example, early responses associated with one region followed by engagement of different regions at later time points. This spatiotemporal progression reflects a dynamic cascade of neural processing that cannot be captured by EEG or fMRI alone.

Building on these early approaches, group-level extensions such as group jICA enable estimation of multimodal components that are consistent across subjects while preserving subject-specific variability. This allows identification of population-level brain networks and facilitates statistical comparisons across groups. As illustrated in Figure 6, group jICA reveals jointly informed components linking fMRI networks with frequency-specific EEG activity. Notably, EEG components projected into source space exhibit spatial patterns that are not present in fMRI-only representations, indicating that fusion methods uncover complementary structure rather than simply combining modalities. This extension highlights a key advance: multimodal fusion not only captures shared dynamics but also improves the spatial interpretability of electrophysiological signals.

More recently, graph neural network models integrating multimodal connectivity identify predictive brain network signatures. As illustrated in Figure 7, multimodal models highlight regions not captured by unimodal analyses, demonstrating that cross-modal integration reveals additional neural features. However, these models introduce challenges in interpretability and data requirements.

However, many classical symmetric fusion methods rely on linear assumptions, whereas the relationship between neural activity and BOLD signals is shaped by complex and nonlinear neurovascular dynamics [3,87]. Some models (notably jICA) also assume shared subject-level mixing across modalities, which, while practical in certain cases, can be restrictive when only part of the underlying neural variance is shared or when one modality contains additional modality-specific structure [123]. Scaling and feature normalization remain important in practice because differences in feature dimensionality (e.g., EEG time–frequency features vs. fMRI voxels/ROIs) can bias the estimation toward one modality if not handled carefully [123]. This constraint can be relaxed by leveraging other models such as parallel ICA.

In contrast, asymmetric fusion uses one modality to constrain analysis in the other. EEG-informed fMRI approaches remain one of the most common asymmetric pipelines, where an EEG-derived regressor (e.g., band power, ERP amplitude, event timing) is used to explain BOLD fluctuations via regression or GLM-based modeling [10,87,90]. These methods often produce highly interpretable outputs, especially when the EEG feature corresponds to a known neurophysiological phenomenon (e.g., alpha suppression, ERP peaks, epileptiform discharges). Similarly, fMRI-informed EEG approaches can constrain EEG source localization or reduce the solution space of the inverse problem, improving spatial interpretability under biophysical assumptions [87,88].

That said, asymmetric fusion is inherently selective: it reveals the aspect of cross-modal coupling targeted by the chosen regressor or constraint but may miss broader shared processes. It can also be sensitive to modeling choices such as HRF shape, lag handling, and event definition, making results more dependent on preprocessing and assumptions than symmetric fusion [87,125]. In practice, symmetric and asymmetric methods should be viewed as complementary: symmetric methods are powerful for discovery and hypothesis generation, while asymmetric methods are often preferred for mechanistic interpretation or clinically actionable outputs.

7.3. Limitations and Open Challenges

Despite progress, several challenges remain persistent. First, interpretation and evaluation remain one of the hardest problems in the field. Unlike unimodal studies where results are directly mapped onto familiar objects (waveforms or activation maps), fusion results often produce latent components or shared representations that require careful interpretation and validation [42,87]. In practice, evaluation is typically multi-criteria: combining cross-modal coupling measures, neurobiological plausibility, and external validation (behavioral or clinical) produces the most convincing evidence. Second, noise and confounds in simultaneous EEG-fMRI continue to shape what fusion methods can recover. Although artifact removal methods have matured, residual gradient and cardioballistic artifacts, head motion, and physiological fluctuations can still bias EEG features and fMRI signals in ways that mimic coupling [3,126]. Third, most fusion frameworks rely on summarized EEG features rather than incorporating the full EEG data. While dimensionality reduction is often necessary for computational tractability, it may discard informative structure. Handling high-dimensional EEG time–frequency–channel data directly within joint decomposition frameworks, without collapsing it into summary measures, remains an open methodological problem. Finally, hemodynamic variability is an unavoidable complication: the BOLD response is delayed and shaped by vascular dynamics, which can vary between regions and individuals, making the direct alignment between EEG events and BOLD responses nontrivial [125,127]. These biophysical mismatches motivate methodological designs that allow flexible coupling, time-lag modeling, or multimodal decompositions that can tolerate partial correspondence [63].

7.4. Future Scope

Scalable group-level fusion frameworks are becoming increasingly important for population-level inference and biomarker discovery [40]. These approaches estimate stable multimodal patterns shared across subjects while retaining subject-specific measures for statistical analysis [10,40,123]. Such designs are particularly useful for comparing clinical cohorts, linking multimodal features to behavioral or symptom measures, and studying variability across individuals.

Another active direction is nonlinear fusion and cross-modal translation using deep learning. Encoder–decoder architectures and generative models can learn complex mappings between EEG and fMRI, potentially allowing EEG to predict spatial fMRI patterns or enabling multimodal latent representations that improve prediction tasks beyond linear models [128].

Flexible coupling models provide an alternative strategy for multimodal integration. Approaches such as soft-coupled tensor decomposition allow EEG and fMRI to share information where the modalities overlap while preserving modality-specific structure when their relationship is partial or nonlinear [59,62,63]. This balance often leads to more robust decompositions and components that remain interpretable across datasets with different acquisition protocols or feature representations.

Although this review focused on data-driven approaches, model-driven and biophysical frameworks remain important for interpreting multimodal findings. Integrating data-driven learning with physiologically informed constraints may help bridge the gap between flexible statistical models and mechanistic understanding of neurovascular coupling [3,87]. Progress along these directions will be important for translating EEG–fMRI fusion from methodological development toward reproducible neuroscience insights and clinically useful biomarkers.

8. Conclusions

Simultaneous EEG-fMRI fusion has evolved from regression-based correspondence toward structured latent-variable modeling capable of identifying shared neural processes across temporal and spatial scales. Data-driven frameworks including ICA-family models, correlation-based methods, tensor decompositions, and emerging nonlinear approaches, have demonstrated that multimodal integration can reveal network-level organization and cross-modal dynamics that are not accessible through unimodal analyses alone. At the same time, challenges remain in interpretation, group-level scalability, modeling spatiotemporal variability, and fully leveraging the richness of high-dimensional EEG data. Although this review focused on data-driven fusion methods applied to simultaneous recordings, other multimodal integration strategies such as model-driven biophysical frameworks, connectivity-based coupling analyses, and studies using non-simultaneous acquisitions have also advanced understanding of EEG-fMRI relationships. Together, these complementary approaches underscore that no single framework is sufficient to capture the full complexity of neurovascular coupling. Continued progress will likely depend on scalable group-based models, flexible coupling strategies that accommodate partial correspondence, and improved evaluation standards. With growing open datasets and methodological rigor, EEG-fMRI fusion is positioned to move beyond methodological development toward reproducible neuroscience insights and clinically meaningful biomarkers.