Optimizing EEG Source Reconstruction with Concurrent fMRI-Derived Spatial Priors.

Reconstructing EEG sources involves a complex pipeline, with the inverse problem being the most challenging. Multiple inversion algorithms are being continuously developed, aiming to tackle the non-uniqueness of this problem, which has been shown to be partially circumvented by including prior information in the inverse models. Despite a few efforts, there are still current and persistent controversies regarding the inversion algorithm of choice and the optimal set of spatial priors to be included in the inversion models. The use of simultaneous EEG-fMRI data is one approach to tackle this problem. The spatial resolution of fMRI makes fMRI derived spatial priors very convenient for EEG reconstruction, however, only task activation maps and resting-state networks (RSNs) have been explored so far, overlooking the recent, but already accepted, notion that brain networks exhibit dynamic functional connectivity fluctuations. The lack of a systematic comparison between different source reconstruction algorithms, considering potentially more brain-informative priors such as fMRI, motivates the search for better reconstruction models. Using simultaneous EEG-fMRI data, here we compared four different inversion algorithms (minimum norm, MN; low resolution electromagnetic tomography, LORETA; empirical Bayes beamformer, EBB; and multiple sparse priors, MSP) under a Bayesian framework (as implemented in SPM), each with three different sets of priors consisting of: (1) those specific to the algorithm; (2) those specific to the algorithm plus fMRI task activation maps and RSNs; and (3) those specific to the algorithm plus fMRI task activation maps and RSNs and network modules of task-related dFC states estimated from the dFC fluctuations. The quality of the reconstructed EEG sources was quantified in terms of model-based metrics, namely the expectation of the posterior probability P(model|data) and variance explained of the inversion models, and the overlap/proportion of brain regions known to be involved in the visual perception tasks that the participants were submitted to, and RSN templates, with/within EEG source components. Model-based metrics suggested that model parsimony is preferred, with the combination MSP and priors specific to this algorithm exhibiting the best performance. However, optimal overlap/proportion values were found using EBB and priors specific to this algorithm and fMRI task activation maps and RSNs or MSP and considering all the priors (algorithm priors, fMRI task activation maps and RSNs and dFC state modules), respectively, indicating that fMRI spatial priors, including dFC state modules, might contain useful information to recover EEG source components reflecting neuronal activity of interest. Our main results show that providing fMRI spatial derived priors that reflect the dynamics of the brain might be useful to map neuronal activity more accurately from EEG-fMRI. Furthermore, this work paves the way towards a more informative selection of the optimal EEG source reconstruction approach, which may be critical in future studies.

Introduction

Electroencephalography (EEG) measures the electrical potential differences between electrodes placed at different scalp sites that are generated by an ensemble of brain cells acting in synchrony. Because of its fairly direct relationship with neuronal activity and its remarkable temporal resolution at the sub-millisecond scale, EEG has proven pivotal for studying both healthy and abnormal human brain function in general, and particularly brain functional connectivity and its dynamics (Niedermeyer and Lopes Da Silva 2005). However, the spatial identification and characterization of the brain networks underlying the electrical potentials measured at the scalp are not possible from these scalp signals alone, given the poor spatial resolution of the EEG at the centimeter scale (Michel et al. 2004). Fortunately, the continuous technological advances of EEG hardware and signal processing techniques now permit a reliable reconstruction of those brain networks (Abreu et al. 2020b), by localizing, taking into account prior information (e.g. fMRI), and estimating the strength of the neural generators responsible for the scalp EEG signals—the so-called EEG source reconstruction (Michel and Murray 2012). However, it is not yet established which EEG source reconstruction (inversion) method to use, nor which prior information is the most useful, considering the application and the context. So, this work is focused on comparing commonly used source reconstruction methods and evaluating the advantages (if any) of adding different types of fMRI information to the reconstruction models.

Reconstructing EEG sources involves a complex pipeline that can be divided into the forward and the inverse problems. The forward problem relates with the estimation of the impact of a given source in the brain on the scalp electrical potentials and is typically solved by building realistic and subject-specific head models from individual structural magnetic resonance images (MRI) using well-defined processing pipelines. By incorporating the 3D localization of the scalp electrodes on these head models, a lead field can then be computed, establishing the relationship between the activity of the different sources in the brain and the signal measured at each electrode (Michel and Brunet 2019). Conversely, the inverse problem relates with determining the sources in the brain that generate a given scalp distribution of electrical potential differences (i.e., EEG topography). Because of the non-uniqueness of its solution, the inverse problem is considered the most challenging, with a plethora of inversion algorithms being available for solving it (Michel et al. 2004). These can be roughly divided into current source density (CSD) estimates and beamformers (Grech et al. 2008; He et al. 2018). Despite the choice between the two types being application-dependent to some extent (Halder et al. 2019), CSD-based algorithms are the most commonly used in the literature; within these, the more recent distributed source localization algorithms are preferable over dipole source localization algorithms, as the latter require prior knowledge regarding the number of sources to be estimated. Several distributed source localization algorithms have been developed, the most common being the minimum norm (MN) solutions (Hämäläinen and Ilmoniemi 1994) and their variations: weighted minimum norm (WMN; De Peralta-Menendez and Gonzalez-Andino, 1998), low resolution electromagnetic tomography (LORETA; Pascual-Marqui et al., 1994), local autoregressive average (LAURA; De Peralta-Menendez et al., 2004), among others (Michel and Brunet 2019). Motivated by the challenging task of defining a ground truth, few studies have dedicated to systematically compare the performance of reconstruction strategies. The prevalent finding is that there is no (universal) optimal solution that encompasses all the contexts and designs, with few studies listing a set of recommendations to guide the choice of the most appropriate method based on the application (Belardinelli et al. 2012; Hedrich et al. 2017; Hincapié et al. 2017; Anzolin et al. 2019; Tait et al. 2021). Particularly, (Hedrich et al. 2017; Hincapié et al. 2017; Tait et al. 2021) have carried a more intensive comparison, however several existing methods were not addressed, e.g. the empirical Bayes beamformer (EBB) and the multiple sparse priors (MSP). Furthermore, these studies were either performed in simulated or resting-state data, and reported inconclusive results (Yao and Dewald 2005; Grova et al. 2006; Bradley et al. 2016; Halder et al. 2019). Similarly, in an attempt to prove the clinical validity of EEG-fMRI, efforts were made, particularly in the context of epilepsy, by using intracranial EEG (icEEG) recordings for validation of EEG source imaging results (Thornton et al. 2010; van Houdt et al. 2013; Vaudano et al. 2013, 2021; Ebrahimzadeh et al. 2021). However, for the chosen reconstruction method their results were only partially validated, i.e., only few patients or regions presented concordant results with both techniques. Although simultaneous icEEG-fMRI can overcome the low sensitivity of scalp EEG and reach good spatial concordance between neuronal electrical and BOLD changes, it comes at the cost of invasive recordings only affordable in specific cases (Vulliemoz et al. 2011; Cunningham et al. 2012; Aghakhani et al. 2015; Chaudhary et al. 2016, 2021; Sharma et al. 2019). Furthermore, these studies did not assess the quality of different source reconstruction algorithms and the effect on concordance measures, which stresses the need to perform systematic comparisons focusing on specific scenarios and considering other less explored methodologies. Moreover, the associated results have not been explicitly validated based on the brain activity of interest; instead, the localization error, spatial spread and percentage of false positives are typically used, as well as the log-evidence and variance explained of the inversion model used when considering Bayesian frameworks (Michel and Brunet 2019).

In addition to choosing the most appropriate inversion method, the choice of a priori information to integrate in these models is also important. This a priori information are assumptions and constraints (priors) that serve to guide the reconstruction and are needed to tackle the non-uniqueness of the inverse problem, which are reflected differently on each inversion algorithm. Their incorporation can be performed using two different approaches: by imposing penalty functions (Valdés-Sosa et al. 2009), or using a Bayesian framework (Trujillo-Barreto et al. 2004; Friston et al. 2008), particularly the parametric empirical Bayesian (PEB; Henson et al., 2010; Phillips et al., 2005). Although less popular, the PEB framework allows to describe a given assumption or constraint explicitly through appropriate postulated prior distributions, which can range from one as in the MN solutions (the identity matrix) to hundreds as in the multiple sparse priors (MSP) algorithm (Friston et al. 2008). This framework is thus extremely flexible for incorporating additional priors obtained from other imaging modalities, which has proved to be beneficial for more efficiently tackling the non-uniqueness of the inverse problem (Liu and He 2006, 2008; Lei et al. 2015). The first studies used brain activation maps obtained from the analyses of task-based functional MRI (fMRI) data (Henson et al. 2010; Lei et al. 2010, 2011, 2012). More recently, the well-known brain networks that emerge from temporally correlated spontaneous fluctuations in the blood-oxygen-level-dependent (BOLD) fMRI signal (the so-called resting-state networks, RSNs) have also been used as spatial priors (Lei 2012). In these studies, the spatial priors were derived from separately acquired fMRI data, which may scale down their potential for guiding the reconstruction of EEG, especially when focusing on spontaneous activity (Abreu et al. 2018). Additionally, task-based and resting-state functional networks are now known to continuously reorganize in response to both internal and external stimuli at multiple time-scales, resulting in temporal fluctuations of their connectivity—the so-called dynamic functional connectivity (dFC) (Hutchison et al. 2013). From dFC fluctuations, a limited, but variable, number of dFC states have been recurrently identified in the literature as the building blocks of brain functional connectivity (dynamics) (Preti et al. 2017), which are hypothesized to be associated with different cognitive, vigilance or pathological brain states (Thompson 2018); however, they have not been considered as potential spatial priors for EEG source reconstruction so far.

Given the increasing relevance of EEG as a brain imaging tool, accurately estimating the underlying brain sources is critical in the study of both healthy and clinical populations. Importantly, no study has so far focused on determining the extent at which the effects of different source reconstruction algorithms and spatial priors differ between groups in clinical studies, with the spatial priors potentially reflecting relevant aspects of the disease under study. This is especially relevant in task-related fMRI studies, which are rapidly increasing in clinical research (Marinazzo et al. 2019). While EEG-fMRI and source imaging reconstruction are highly applied in the context of epilepsy (Gotman and Pittau 2011; Lei et al. 2015; van Graan et al. 2015), applying these techniques to investigate other neurological and psychiatric diseases in which altered connectivity is suspected can potentially result in highly useful clinical applications. Particularly, Multiple sclerosis (MS) is a disconnection disease that is due to structural damage but also functional connectivity alterations, with both EEG (with its high temporal resolution) and fMRI (with its high spatial resolution) representing gold-standard techniques to investigate it (Gschwind et al. 2016; Tahedl et al. 2018). By leveraging the high temporal resolution from EEG and spatial resolution from fMRI, the underlying temporally- and spatially-resolved EEG sources with fMRI-derived information will provide robust connectivity measures that might help to understand the pathophysiology of the disease and serve as a tool for reliable assessment of disease progression. Moreover, while it is very common to explore connectivity measures in patients with MS during resting state (Sbardella et al. 2015), task-designs target brain regions and networks that show distinct properties than in resting-state (Di et al. 2013). Thus, task designs may have a crucial role in describing the functioning of the brain, in highlighting specific connectivity changes, and thus in understanding this disease better.

Considering the present limitations described in this section, here we compared four different inversion algorithms (MN, LORETA, empirical Bayes beamformer, EBB; and MSP), each with two different sets of additional fMRI-derived spatial priors (activation maps and RSNs, with and without including dFC states) on EEG data collected concurrently with fMRI at 3 T from 6 multiple sclerosis (MS) patients and 7 healthy subjects performing visual perception tasks and during rest. The quality of the reconstructions was quantified through the expectation of the posterior probability P(model|data), obtained from the log-evidence, and variance explained of the associated models, and in terms of the overlap between EEG source components and brain regions of interest associated with the tasks and RSN templates.

Materials and Methods

Participants

Six MS patients (mean age: 30 ± 8 years; 2 males) and seven demographically matched healthy subjects (mean age: 30 ± 6 years; 3 males) were recruited. The patients were selected by the clinical team at the Neurology Department of the University Hospital of Coimbra, and met the criteria for MS diagnosis according to McDonald Criteria (Thompson et al. 2018). All participants had normal or corrected-to-normal vision.

Experimental Protocol

The imaging session was performed at the Portuguese Brain Imaging Network (Coimbra, Portugal) and consisted of four functional runs: first, a functional localizer of the human middle temporal area (hMT + /V5, a low level visual area well-known to respond to simple motion patterns), followed by two runs of biological motion (BM) perception, and one final resting-state run.

The localizer run consisted of 10 blocks of 18 s, with each block comprising three periods: the first was a fixation period marked by a red cross positioned at the center of the screen for 6 s. During the second period, a 6 s pattern of stationary dots was shown, followed by the third (and final) period during which the dots were moving towards and away from a central fixation cross at a constant speed (5 deg/sec) for 6 s.

Biological motion stimuli were built based on human motion capture data collected at 60 Hz, comprising 12 point-lights placed at the main joints of a male walker. Each BM perception run consisted of 12 blocks of 40 s: 4 or 5 blocks (depending on the starting block) of the point-light walker facing rightwards or leftwards (body blocks), 4 or 5 blocks showing only the point-light located at the right ankle and moving rightwards of leftwards (foot blocks), and 3 blocks of the original 12 point-lights randomly positioned across the y axis, while maintaining their true trajectory across the x axis (random blocks). A total of 9 body, 9 foot and 6 random blocks were then collected during the two BM perception runs.

During the resting-state run, the participants were instructed to relax and only fixate a red cross positioned at the center of the screen.

EEG-fMRI Data Acquisition

Imaging was performed on a 3 T Siemens MAGNETOM Prisma Fit MRI scanner (Siemens, Erlangen) using a 64-channel RF receive coil. In order to minimize head motion and scanner noise related discomfort, foam cushions and earplugs were used, respectively. The functional images were acquired using a 2D simultaneous multi-slice (SMS) gradient-echo echo-planar imaging (GE-EPI) sequence (6 × SMS and 2 × in-plane GRAPPA accelerations), with the following parameters: TR/TE = 1000/37 ms, voxel size = 2.0 × 2.0 × 2.0 mm3, 72 axial slices (whole-brain coverage), FOV = 200 × 200 mm2, FA = 68°, and phase encoding in the anterior–posterior direction. The start of each trial was synchronized with the acquisition of the functional images. A short EPI acquisition (10 volumes) with reversed phase encoding direction (posterior-anterior) was also performed prior to each fMRI run, for image distortion correction. Whole-brain, 1 mm isotropic structural images were acquired using a T1-weighted 3D gradient-echo MP2RAGE sequence.

The EEG signal was recorded using the MR-compatible 64-channel NeuroScan SynAmps2 system and the Maglink™ software, with a cap containing 64 Ag/AgCl non-magnetic electrodes positioned according to the 10/10 coordinate system, a dedicated electrode for referencing placed close to the Cz position, and two electrodes placed on the back for electrocardiogram (EKG) recording. Electrode impedances were kept below 25 kΩ. EEG, EKG and fMRI data were acquired simultaneously in a continuous way, and synchronized by means of a Syncbox (NordicNeuroLab, USA) device. EEG and EKG signals were recorded at a sampling rate of 10 kHz, synchronized with the scanner’s 10 MHz clock. No filters were applied during the recordings. The helium cooling system was not turned off, as it may carry the associated risk of helium boil-off in certain systems (Mullinger et al. 2008), and thus is not permitted in some clinical sites as the one used in this study. Respiratory traces were recorded at 50 Hz with a respiratory cushion from the physiological monitoring unit of the MRI system.

For each participant, 197 fMRI volumes were acquired during the localizer run, yielding 3.20 min of duration. The two BM runs had approximately 8.37 min, thus comprising 507 volumes each. The final resting-state run had approximately 8.08 min, corresponding to 485 volumes.

MRI Data Analysis

The main steps of the processing pipeline for deriving fMRI spatial priors (described here) and subsequently use them in EEG source reconstruction (described in the next section), as well as the metrics proposed for quantifying the quality of the source reconstruction, are depicted in Fig. 1.

Fig. 1

Schematic diagram of the processing pipeline. The pre-processed fMRI data is submitted to three different analyses in order to derive three types of fMRI spatial priors for EEG source reconstruction: (1) identification of RSNs through spatial ICA; (2) mapping of the task-related activity through GLM; and (3) by estimating the dFC fluctuations with phase coherence and the associated dFC states with dictionary learning, dFC state modules were obtained using the Louvain modularity algorithm. The covariance components (CCs) associated with these spatial priors were then included in several inversion algorithms, whose reconstruction quality was assessed by the expectation of the posterior probability P(model|data) and variance explained of the associated models, and by the overlap of EEG source components (obtained through spatial ICA applied to the source reconstructed EEG) with ROIs and RSN templates

Pre-processing

The first 10 s of data were discarded to allow the signal to reach steady-state. Subsequently, slice timing and motion correction were performed using FSL tool MCFLIRT (Jenkinson et al. 2002), followed by a B0-unwarping step with FSL tool TOPUP (Andersson et al. 2003) using the reversed-phase encoding acquisition, to reduce EPI distortions. The distortion-corrected images were then corrected for the bias field using FSL tool FAST (Zhang et al. 2001), and non-brain tissue was removed using FSL tool BET (Smith 2002). Nuisance fluctuations (including physiological noise) were then removed by linear regression using the following regressors (Abreu et al. 2017): (1) quasi-periodic BOLD fluctuations related to cardiac and respiratory cycles were modeled by a fourth order Fourier series using RETROICOR (Glover et al. 2000); (2) aperiodic BOLD fluctuations associated with changes in the heart rate as well as in the depth and rate of respiration were modeled by convolution with the respective impulse response functions (as described in Chang et al. 2009); (3) the average BOLD fluctuations measured in white matter (WM) and cerebrospinal fluid (CSF) masks (obtained as described below); (4) the six motion parameters (MPs) estimated by MCFLIRT; and (5) scan nulling regressors (motion scrubbing) associated with volumes acquired during periods of large head motion; these were determined using the FSL utility fsl_motion_outliers, whereby the DVARS metric proposed in (Power et al. 2012) is first computed, and then thresholded at the 75th percentile plus 1.5 times the inter-quartile range. Finally, a high-pass temporal filtering with a cut-off period of 100 s was applied, and spatial smoothing using a Gaussian kernel with full width at half-maximum (FWHM) of 3 mm was performed.

For each subject, WM and CSF masks were obtained from the respective T1-weighted structural image by segmentation into gray matter, WM and CSF using FSL tool FAST (Zhang et al. 2001). The functional images were co-registered with the respective T1-weighted structural images using FSL tool FLIRT, and subsequently with the Montreal Neurological Institute (MNI) (Collins et al. 1994) template, using FSL tool FNIRT (Jenkinson and Smith 2001; Jenkinson et al. 2002). Both WM and CSF masks were transformed into functional space and were then eroded using a 3 mm spherical kernel in order to minimize partial volume effects (Jo et al. 2010). Additionally, the eroded CSF mask was intersected with a mask of the large ventricles from the MNI space, following the rationale described in (Chang and Glover 2009).

Each participant’s structural image was parceled into non-overlapping regions of the cerebrum according to the automated anatomical labeling (AAL) atlas (Tzourio-Mazoyer et al. 2002). These parcels were co-registered to the participant’s functional space, and the pre-processed BOLD data were then averaged within each parcel.

fMRI Priors for EEG Source Reconstruction

From the pre-processed fMRI data, several potential priors for EEG source reconstruction were subsequently extracted (procedures described next), namely: resting-state networks for all runs, and task-related activity maps and dynamic functional connectivity (dFC) states for the task runs only.

Identification of Resting-State Networks

The pre-processed fMRI data were submitted to a group-level probabilistic spatial ICA (sICA) decomposition using the FSL tool MELODIC (Beckmann and Smith 2004), whereby the data of each run for all participants is temporally concatenated prior to the sICA step, as recommended in the MELODIC’s guide for the identification of RSNs (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/MELODIC). The optimal number of independent components (ICs) was automatically estimated based on the eigenspectrum of its covariance matrix (Beckmann and Smith 2004), with an average of approximately 40 ICs across runs.

An automatic procedure for the identification of well-known RSNs was then applied, in which the spatial maps of the ICs (thresholded at Z = 3.0) were compared with those of the 10 RSN templates described in (Smith et al. 2009), in terms of spatial overlap computed as the Dice coefficient (Dice 1945). For each template, the IC map yielding the highest Dice coefficient was determined as the corresponding RSN. In the cases of non-mutually exclusive assignments, the optimal assignment was determined by randomizing the order of the RSN templates (a maximum of 10,000 possible combinations were considered, for computational purposes), and then sequentially, and mutually exclusively, assigning them to the IC maps based on their Dice coefficient. The assignment with the highest average Dice coefficient across all RSN templates was then deemed optimal, yielding the final set of RSNs: three visual networks (RSN 1–3), the default mode network (DMN, RSN4), a cerebellum network (RSN5), a motor network (RSN6), an auditory network (RSN7), the salience network (RSN8), a right language network (RSN9) and a left language network (RSN10).

These 10 subject- and run-specific RSNs were then used as spatial priors for the reconstruction of sources of EEG collected on all four runs. RSNs were considered in the task runs because these have been shown to be also present in task-based studies (Di et al. 2013; Cole et al. 2016).

hMT + and BM-related Activity Mapping

For the purpose of mapping hMT + /V5 from the localizer run, and the regions involved in the BM perception task from the other two runs, a general linear model (GLM) framework was used. For the localizer, two regressors representing the periods showing dots (stationary and moving) were considered. These regressors were built based on unit boxcar functions with ones during the respective periods, and zeros elsewhere. Similarly, three regressors representing the body, foot and random blocks of the BM runs were built for analyzing the BM runs, with the regressors also based on unit boxcar functions. All regressors were convolved with a canonical, double-gamma hemodynamic response function (HRF). The run-specific, HRF-convolved regressors were then included in a GLM that was subsequently fitted to the associated fMRI data using FSL tool FILM (Woolrich et al. 2001). The hMT + /V5 regions were identified from the localizer run by contrasting the moving and the stationary dots periods, whereas the areas associated with BM perception were mapped according to the following contrasts: body—random, foot—random, and body—foot. Voxels exhibiting significant changes within these contrasts were identified by cluster thresholding (voxel Z > 2.5, cluster p < 0.05).

In this way, a single spatial prior is obtained for reconstructing the sources of the EEG collected during the localizer run. In contrast, three spatial priors (one for each contrast) are made available for each of the two BM runs.

Dynamic Functional Connectivity Analysis

The dFC analysis described here was only performed on the fMRI data collected during the task runs (localizer and two BM runs), since its purpose was to objectively identify a small set of dFC states associated with the tasks, and to use them as spatial priors in the subsequent reconstruction of EEG sources. The dFC was estimated using the phase coherence (PC) method, which allows to compute the dFC for each fMRI sample; the loss in temporal resolution and the ambiguous selection of a window size, both inevitable in conventional sliding-window correlation approaches (Preti et al. 2017), are thus avoided (Glerean et al. 2012). For the PC method, a second-order Butterworth band-pass filter in the range of 0.01–0.1 Hz was first applied to the parcel-averaged BOLD signals. The instantaneous phase, θ, of the filtered signal, n, was then estimated using the Hilbert transform (Cabral et al. 2017; Figueroa et al. 2019). One of the advantages of this method is to obtain a connectivity matrix for each TR (unlike Pearson correlation, which requires a sliding window). Also, the real part of the signals does not depend on their amplitude, and therefore it is a measure less sensitive to the presence of noise, namely movement that creates large amplitude variations in the fMRI signal. For each participant, the dFC matrix (N = 90 brain regions from the AAL atlas, and T is the number of fMRI samples, which depends on the run under analysis) was computed for each pair of parcels, n and p, and at each fMRI sample t, using the equation:

For each run and participant individually, the matrix C was then submitted to the leading eigenvector dynamics analysis (LEiDA) (Cabral et al. 2017; Figueroa et al. 2019; Lord et al. 2019), with the purpose of reducing the dimension of each temporal entry of C (N × N) by only considering the associated leading eigenvector (of dimension N), while nonetheless explaining most of the variance (> 50% in all cases, and up to 90%) (Lord et al. 2019). This step yielded the reduced dFC matrix , with the columns (t = 1, …, T) representing the leading eigenvectors, and the rows indicating the parcels. Each eigenvector is composed by elements with positive and/or negative signs; if all positive, a global mode is governing the parcel-averaged BOLD signals where all the associated phases point in the same direction with respect to the orientation defined by the eigenvector (Figueroa et al. 2019). If the elements of the eigenvector have different signs, the parcels can be grouped into two networks according to their sign (positive or negative) in the eigenvector. The magnitude of the elements indicates how strongly the associated parcel belongs to its assigned network (Newman 2006).

For the identification of dFC states, an l1-norm regularized dictionary learning (DL) approach was employed, following the methodology proposed in (Abreu et al. 2019). Briefly, this can be formulated as the matrix factorization problem , where and represent the dFC states and associated weight time-courses (i.e. contribution of each dFC state to reconstruct C at each time point), respectively; and k is the number of dFC states. These are estimated by solving the optimization problem given by:so that the reconstruction error of C, , is minimized; denotes the Frobenius norm of a matrix. The estimation of D and A was performed using the algorithms implemented in the MATLAB® toolbox SPArse Modeling Software (SPAMS, Mairal et al., 2010). The sparsity of the solutions was controlled by a non-negative parameter λ on an l1-norm regularization framework. The number of dFC states k was varied between from 5 to 10 in unit steps, and λ between ten values from 1 to 0.1259 in decreasing exponential steps.

The optimal k and λ values were jointly determined with the dFC states to be considered as spatial priors in the EEG source reconstruction. This was achieved by first computing the Pearson correlation between the contrasts defined for identifying the activation maps (one for the localizer run, and three for the BM runs) and the dFC weight time-courses in A, for all possible combinations of k and λ. For the localizer run, the dFC state exhibiting the highest correlation across dFC states, and combinations of k and λ, was deemed as task-related. For the BM runs, the dFC state exhibiting the highest correlation across contrasts, dFC states and combinations of k and λ was first identified. For the optimal combination of k and λ, the most correlated dFC states associated with the remaining contrasts were then determined. In cases where multiple contrasts were associated with the same dFC state, only that state was considered for the subsequent analyses.

The dFC states of interest were then finally submitted to a modularity analysis, with the purpose of identifying modular (or community) structure in those states. Because the dFC states are column vectors, rather than square matrices representing a connectivity matrix as required for the modularity analysis, such connectivity matrix of a given dFC state was first reconstructed by computing the outer product of the dFC state, (Cabral et al. 2017). The Louvain algorithm as implemented in the Brain Connectivity Toolbox was then applied to the reconstructed connectivity matrices of the dFC states of interest (Rubinov and Sporns 2010). This algorithm considers both the positive and negative weights of the unthresholded connectivity matrix, thus avoiding the ambiguous selection of a threshold as required in conventional modularity algorithms. Each of the N parcels is then assigned a label indicating which module the parcel belongs to. The network modules were then projected into binary 3D spatial maps to be used as spatial priors in the EEG source reconstruction, by identifying the voxels belonging to parcels (according to the AAL atlas used for parceling the brain) assigned to the same modules. The number of modules automatically identified by the Louvain algorithm was between 2 and 4 in all cases; in this way, the number of spatial priors built from this analysis varied according to the number of contrasts (run) and modules identified, with a maximum of 1 [contrast] × 4 [modules] = 4 for the localizer run, and 3 [contrasts] × 4 [modules] = 12 for the BM runs.

EEG Data Analysis

EEG data underwent gradient artifact correction on a volume-wise basis using a standard artifact template subtraction (AAS) approach (Allen et al. 2000) using the FMRIB tools implemented as a plug-in of the EEGLAB toolbox (Delorme and Makeig 2004). The pulse artifact was removed using the method presented in (Abreu et al. 2016), whereby the EEG data is first decomposed using independent component analysis (ICA), followed by AAS to remove the artifact occurrences from the independent components (ICs) associated with the artifact. The corrected EEG data is then obtained by reconstructing the signal using the artifact-corrected ICs together with the original non-artifact-related ICs.

After the removal of the MR-induced artifacts, EEG data was then submitted to some of the routines of the automatic pipeline (APP) for EEG pre-processing described in (da Cruz et al. 2018), namely: (1) re-referencing to a robust estimate of the mean of all channels; (2) removal and interpolation of bad channels; and (3) removal of bad epochs of 1 s (matching the TR of the fMRI data). An additional ICA step was then performed with the purpose of removing additional sources of EEG artifacts; these were identified using the ICLabel algorithm (Pion-Tonachini et al. 2019), implemented as a plug-in of the EEGLAB toolbox (Delorme and Makeig 2004). The classification provided by ICLabel is based on a previously trained model with a large EEG dataset collected outside the MR scanner, rendering this algorithm sub-optimal for our dataset. To cope with this, all ICs were visually inspected in order to validate, and eventually correct (for both false positives and negatives), the classification outputs of ICLabel. Finally, the EEG data was down-sampled to 500 Hz and band-pass filtered to 1–30 Hz.

Source Reconstruction

The pre-processed EEG data from all runs was then submitted to several EEG source reconstruction procedures implemented in SPM12 (https://www.fil.ion.ucl.ac.uk/spm/). To reduce the computational load, the EEG data was further downsampled to a sampling rate of 60 Hz (two times the highest frequency component of the data) using a polyphase anti-aliasing FIR filter as implemented in EEGLAB (Delorme and Makeig 2004).

The Forward Problem

A realistic head model was built by first segmenting each participant’s structural image into 3 tissue labels (brain, scalp and skull), and computing the deformation field needed to co-register the structural images into an MNI template. The individual meshes were then obtained by applying the inverse of this deformation field to the canonical meshes derived from the MNI template; meshes with 8196 vertices were considered. The electrode positions were co-registered to the scalp compartment by first considering their standard positions (in the 10/10 coordinate system), and then manually adjusting them to match the distortions clearly observed on the structural images. A realistically shaped volume conduction model was estimated using a boundary element model (BEM) with three layers (scalp, inner skull and outer skull). 8196 source dipoles were placed at the vertices of a cortical surface also derived from the MNI template and transformed into the structural image. The leadfield matrix was then estimated, mapping each possible dipole configuration onto a scalp potential distribution.

The Inverse Problem

The inverse problem was solved using a Parametric Empirical Bayesian (PEB) framework as implemented in SPM12, which can be formulated as (López et al. 2014):where is the EEG data with C channels (64 in this case) and T time samples (depends on the run under analysis); is the leadfield matrix (D is the number of dipoles, 8196 in this case); and is the unknown source dynamics at each dipole. represents the multivariate Gaussian probability distribution, and T the temporal correlations (known and fixed). The terms and denote the noise at the channel and source spaces, with covariance matrices and , respectively. Channel noise is typically assumed to be uniform across channels, and therefore can be defined as , with hC the channel noise variance and the identity matrix. The source space covariance matrix CD assumes the form:where represents different types of covariance components (CCs) reflecting prior knowledge on the sources to be reconstructed, and the unknown hyperparameter denoting its relative importance. These hyperparameters work as regularization parameters in ill-posed problems such as the EEG inverse problem, and were estimated using a restricted maximum likelihood (ReML) algorithm that uses as cost function the log-evidence of the model. Commonly used source inversion algorithms can then be derived from Eq. 3 by defining the CCs that appropriately reflect their assumptions. For instance, MN solutions assume that all dipoles have the same variance and no covariance; therefore, only one CC is defined as , with hD the source noise variance and the identity matrix.

In this work, we tested four source inversion algorithms: MN, LORETA, EBB and MSP; their derivations from Eq. 3 and associated CCs are thoroughly presented in (López et al. 2014). Additionally to those specific to a given algorithm, other CCs estimated from the fMRI-derived spatial priors (RSNs, activation maps and task-based dFC states) were considered (the procedures for their estimation are briefly described below). Specifically, all four inversion algorithms were tested using three different sets of CCs: the simplest set S1, with only CCs specific to the algorithm: (2) a larger set S2 comprising S1 and CCs from RSNs and activation maps (the latter for task runs only); and (3) the largest set S3 comprising S2 and CCs from the modules of the task-related dFC states (hence only tested on EEG collected from task runs). A total of 4 [inversion algorithms] × 3 [sets of CCs] = 12 reconstructions of EEG sources S were then performed for each subject and run (only 8 for the resting-state runs).

Estimation of Covariance Components from fMRI-Derived Spatial Priors

CCs were estimated from the fMRI-derived spatial priors by first transforming them into binary priors. These 3D binary spatial priors were then projected onto the 2D cortical surface using nearest-neighbor interpolation (Henson et al. 2010), and smoothed using the Green’s function G of the cortical mesh adjacency matrix , (Harrison et al. 2007). The entries of M, mij, are 1 if vertices i and j of the cortical surface are neighbors (within a defined radius) and 0 otherwise; here, a radius of 8 vertices and a smoothing parameter of were selected according to (Friston et al. 2008). The CCs of the smoothed (and projected) spatial priors are then obtained by computing their covariance matrices, i.e., their outer product. These procedures are illustrated in Fig. 2.

Fig. 2

Deriving covariance components (CCs) from fMRI spatial priors. The 3D fMRI spatial priors are first binarized, projected onto the 2D cortical surface using nearest-neighbor interpolation and smoothed using the Green’s function. The associated CCs are then obtained by computing the outter product. For visualization purposes, the temporally reduced CCs are illustrated, by applying the same temporal projector considered when reducing the EEG data prior to its reconstruction

Source Reconstruction Quality

EEG Source Components

Following the rationale of previous studies (Liu et al. 2017, 2018; Abreu et al. 2020b), a spatial ICA step similar to that applied to the fMRI data for identifying RSNs was then performed on the reconstructed source dynamics S, with the purpose of separating those potentially associated with RSNs and/or other regions of interest in our tasks. This can be formulated as:where is the mixing matrix, with each column the time-course of the source component (SC) i; and represents the spatial maps in the source space associated to each of the I SCs. Because the EEG data is submitted to a temporal reduction step prior to solving the inverse problem in order to reduce noise while guaranteeing a temporally continuous estimation of sources (López et al. 2014), the rank of S is reduced accordingly, being then defined an upper bound on the number of SCs to be estimated. Such maximum allowed number of SCs was then estimated, which was between 50 and 60 in all cases. Finally, the SCs were converted into z-scores, and the deformation field estimated while solving the forward problem was applied to transform them from the source space into the MNI space.

Quality Metrics

Model-based quality metrics were first considered, namely the variance explained (VE) of the reconstructed EEG data relative to the actual EEG data Y (see Eq. 3), and the expectation of the posterior probability P(model | data) of the inversion models. The latter is obtained from the associated log-evidence values of all models, for all subjects, as described in (Rigoux et al. 2014), and reflects the probability of obtaining a given model when randomly selecting a subject. These probability values were normalized to sum to one, over the models under analysis. Other quality metrics reflecting more directly the presence of neuronal activity of interest in the SCs were also considered, as described next.

First, because the perception of motion in general, and of biological motion in particular, is known to elicit certain brain regions (Chang et al. 2018), the following four spherical regions of interest (ROIs) of 10 mm centered at specific MNI coordinates (indicated in square brackets) were considered: anterior insula (aINS) at [± 36, 24, 2], extrastriate body area (EBA) at [left –46, –75, –4; right 47, –71, –4], fusiform body area (FBA) at [left –38, –38, –27; right 43, –43, –28], and fusiform gyrus (FFG) at [± 42, –56, –14]. Four additional task-related brain regions (Chang et al. 2018) were obtained from FSL atlases (threshold applied to the probability maps is indicated in square brackets), namely: inferior frontal gyrus (IFG) [0.25], posterior superior temporal sulcus (pSTS) [0.25], visual area V3 [0.25] and visual area hMT + /V5 [0.10]. After binarizing the ROIs and the SC maps, the Dice coefficient d, and the proportion of the ROIs contained in the SC maps p, were then quantified according to (Dice 1945):where N and N denote the number of non-zero voxels in the ROIs and SC maps, respectively, and N the number of overlapping non-zero voxels between the two images; both measures range from 0 (no overlap) to 1. These same two measures, d and p, were also computed between the SC maps and 10 RSN templates described in (Smith et al. 2009), in order to assess which, if any, SCs represented RSNs (similar to the identification of RSNs on fMRI data described previously).

All these measures were computed for each subject, run, inversion algorithm, set of covariance components (CCs), SC maps and maps of interest (8 ROIs and 10 RSN templates). Because only a subset of the SC maps is expected to be associated with those maps of interest, the SC map yielding the highest dice coefficient for each map was identified, and the associated and maximum values kept for subsequent analyses. The and values were further summarized by computing their average within each map type (ROIs and RSN templates), thus yielding the final set of 13 [subjects] × 4 [runs] × 4 [inversions] × 3 [sets of CCs] × 2 [map types] = 1440 values of and .

Statistical Analysis

The main effects of the population group (MS patients and healthy subjects), inversion algorithm, the set of CCs and the type of map of interest, as well as interaction effects, were evaluated by means of a 4-way repeated measures Analysis of Variance (ANOVA) for the VE, and measures treated separately as the dependent variables. Multiple comparisons between the inversion algorithms, sets of CCs and interactions between the two were performed by means of a post-hoc statistical test with the Tukey–Kramer correction. A level of statistical significance p < 0.05 was considered.

Results

In this work, the quality of EEG source reconstruction provided by the different combinations of (four) inversion algorithms and (three) sets of CCs, was first evaluated in terms of the posterior P(model|data) and VE of the associated models, which are commonly considered in PEB frameworks. Because no significant differences were observed between population groups (healthy subjects and MS patients) the VE values shown in Table 1 were averaged across participants; the values associated with the three visual perception task runs (hMT + /V5 functional localizer and two BM runs) were also averaged. The combination MSP + S1 (with S1 containing only CCs specifically associated with the inversion algorithm) yielded the highest P(model|data) and the highest VE, only followed by LORETA + S2 and LORETA + S3. The ANOVA of the VE values revealed significant main effects of the inversion algorithms and sets of CCs, as well as a significant interaction. The post-hoc test showed no statistically significant differences between inversion algorithms, and the set S1 was significantly better than the sets S2 and S3; as expected, the combination MSP + S1 performed significantly better than other five combinations using the MN or EBB as inversion algorithms and sets S2 or S3.

Table 1

Posterior P(model|data) and average VE values across participants, and across three visual perception task runs, for all combinations of inversion algorithms and sets of covariance components

Sets of CCs	Inversion algorithms	Task runs (Localizer + BMs)		Resting-state runs
		P [%]	VE [%] (± std)	P [%]	VE [%] (± std)
S1	MN	4.0	79.7 ± 14.6	4.8	84.1 ± 11.2
	LORETA	4.0	79.5 ± 14.5	4.8	83.9 ± 11.3
	EBB	4.0	79.4 ± 14.6	4.8	83.2 ± 11.8
	MSP	56.0	84.5 ± 10.5	66.4	87.6 ± 8.3
S2	MN	4.0	81.8 ± 13.7	4.8	85.6 ± 10.6
	LORETA	4.0	81.9 ± 13.7	4.8	85.5 ± 10.6
	EBB	4.0	79.4 ± 14.7	4.8	83.2 ± 11.8
	MSP	4.0	63.9 ± 20.4	4.8	74.9 ± 19.5
S3	MN	4.0	80.7 ± 13.7	NA	NA
	LORETA	4.0	80.7 ± 13.7	NA	NA
	EBB	4.0	78.2 ± 14.8	NA	NA
	MSP	4.0	60.3 ± 15.5	NA	NA

Values in bold represent the best across inversion algorithms for each CC set, and values in bold italics represent the overall best (across inversion algorithms and CC sets)

In order to directly reflect the presence of neuronal activity of interest in the EEG source components (SCs), the reconstruction quality was then quantified in terms of the overlap of SCs with the 8 ROIs and the 10 RSN templates from (Smith et al. 2009). This is illustrated in Fig. 3, showing a considerable overlap (in terms of and ) of two SCs with the EBA mask and the visual RSN1 template, for the first BM run of a given healthy subject. Consistently with the P(model|data) and VE values, the and values were not statistically significantly different between population groups, and thus were averaged across participants and across task runs; these are depicted in Table 2. When considering only the CCs specific to the inversion algorithms (set S1), EBB yields the best results in all cases, and is the overall best (across sets of CCs) in terms of for the resting-state run. However, by combining S1 with RSNs and activation maps (set S2), MN achieves the highest and values for both types of runs, and the overall highest values (across sets of CCs) in terms of for the resting-state run. For the task runs, the largest set of CCs including the dFC state modules (set S3) exhibits the overall best source reconstruction. Similarly to the statistical analysis of VE, the ANOVA of the and values revealed significant main effects of the inversion algorithms and sets of CCs, as well as a significant interaction. For the values, the post-hoc statistical tests showed that MN and EBB inversion algorithms performed significantly better than LORETA and MSP, and that using the sets S2 or S3 was significantly better than only considering the set S1. The latter observation was also true for the values, although in this case it was the MN and MSP inversion algorithms that yielded significantly better results than LORETA and EBB. The combinations EBB + S2 and MSP + S3 exhibited significantly higher and values, respectively, than a subset of combinations including the MN and LORETA algorithms, and the set S1.

Fig. 3

Illustration of the overlap between two EEG SCs (in red-yellow) and A the EBA mask (in blue) and B a visual RSN (in blue-light blue) from (Smith et al. 2009). The dice coefficient d and the proportion of the ROIs contained in the respective SCs are also depicted

Table 2

Average and values across participants, and across three visual perception task runs, for all combinations of inversion algorithms and sets of covariance components

Sets of CCs	Inversion algorithms	Task runs (Localizer + BMs)		Resting-state runs
		d* (± std)	pRS∗(± std)	d* (± std)	pRS∗(± std)
S1	MN	0.15 ± 0.05	0.22 ± 0.10	0.15 ± 0.05	0.22 ± 0.10
	LORETA	0.12 ± 0.05	0.14 ± 0.07	0.12 ± 0.04	0.15 ± 0.06
	EBB	0.15 ± 0.05	0.24 ± 0.11	0.16 ± 0.06	0.25 ± 0.11
	MSP	0.14 ± 0.05	0.18 ± 0.09	0.12 ± 0.04	0.16 ± 0.09
S2	MN	0.15 ± 0.05	0.24 ± 0.12	0.16 ± 0.05	0.24 ± 0.11
	LORETA	0.15 ± 0.05	0.23 ± 0.12	0.15 ± 0.04	0.22 ± 0.12
	EBB	0.15 ± 0.05	0.21 ± 0.10	0.15 ± 0.05	0.23 ± 0.11
	MSP	0.14 ± 0.05	0.23 ± 0.11	0.14 ± 0.05	0.22 ± 0.09
S3	MN	0.15 ± 0.05	0.24 ± 0.12	NA	NA
	LORETA	0.16 ± 0.05	0.24 ± 0.12	NA	NA
	EBB	0.15 ± 0.05	0.23 ± 0.11	NA	NA
	MSP	0.12 ± 0.07	0.35 ± 0.10	NA	NA

Values in bold represent the best across inversion algorithms for each CC set, and values in bold italics represent the overall best (across inversion algorithms and CC sets)

Discussion

In this work, we aimed at optimizing the reconstruction of EEG sources by considering spatial priors derived from concurrently acquired fMRI data when solving the inverse problem, coupled with a systematic comparison of different inversion algorithms and sets of covariance components (CCs) reflecting those spatial priors, on a parametric empirical Bayesian (PEB) framework. The quality of the source reconstructions was quantified in terms of PEB-based metrics (the expectation of the posterior probability, P(model|data); and the variance explained of the respective inversion models, VE), and physiologically-based metrics (the overlap of EEG source components with ROIs and RSN templates representative of brain activity of interest), the latter directly reflecting the presence of neuronal activity.

EEG Source Reconstruction Quality

Under a PEB framework, four inversion algorithms were tested here (MN, LORETA, EBB and MSP) for reconstructing sources from real EEG data collected from participants performing visual perception tasks and during rest and considering three different sets of CCs. We found that depending on the type of quality metric (PEB- or physiologically-based), different conclusions could be taken, in line with the notion that these metrics carry qualitatively distinct information. In terms of the PEB-based metrics (P(model | data) and VE), using the set consisting only of CCs specific to the inversion algorithms (S1) always yielded significantly better results than the CC sets including fMRI spatial priors (S2 and S3, comprising activation maps and RSNs, with or without dFC state modules, respectively). In contrast, by considering S2 and S3, the overlap of EEG source components (SCs) with the ROIs and RSN templates (measured by the dice coefficient d, and the proportion of the ROIs/RSNs contained in the SC maps, p) significantly surpassed that of S1. On the one hand, these contrasting results evidence the underlying optimization procedure used here, combining the PEB framework with the restricted maximum likelihood (ReML) algorithm for estimating the hyperparameters associated with each CC (Phillips et al. 2005; Henson et al. 2010). In fact, adding fMRI spatial priors drastically increases model complexity, which is penalized by ReML, and thus may explain the best PEB-based metrics when considering the more parsimonious inversion models (López et al. 2014). On the other hand, assessing the source reconstruction quality with metrics reflecting more directly the presence of neuronal activity of interest revealed that the information contained on the fMRI spatial priors is pivotal, suggesting that increasing model complexity in this way is needed for EEG SCs to contain such activity of interest, which is of the utmost interest for any subsequent analyses. Accordingly, the usefulness of fMRI spatial priors on EEG source reconstruction has already been shown in previous studies, with the addition of task-based activation maps (Henson et al. 2010; Lei et al. 2010, 2011, 2012) or RSNs (Lei 2012) similar to those considered here improving the reconstructions. These however have only been compared in terms of conventional quality metrics, without taking explicitly into account the neuronal activity of interest. These observations highlight the relevance of using multiple quality metrics expressing different aspects of the reconstructed sources for more appropriately characterizing them, and consequently better informing and improving the selection of the reconstruction approach.

Regarding the optimal inversion algorithm, we found that no statistical differences were observed when comparing the P(model|data) and VE values, whereas MN and EBB yielded significantly higher d values than LORETA and MSP; MN and MSP achieved significantly higher p values than EBB and LORETA. In contrast with the remaining inversion algorithms, LORETA is known for its low-resolution solutions (Michel and Murray 2012; Michel and Brunet 2019), which may render it inappropriate for localizing sources specifically associated with the limited number of rather small brain regions known to be involved in the tasks used in this study, and thus explaining its poorest performance (Halder et al. 2019). Concordant observations have been reported on previous comparison studies on simulated data (Yao and Dewald 2005; Grova et al. 2006; Bradley et al. 2016; Halder et al. 2019), although no differences in performance were shown between MN and LORETA on real magnetoencephalography (MEG) or high-density EEG data (Hedrich et al. 2017). Similarly to MN, MSP has also been shown to provide solutions with high resolution (measured by the focal activation, for instance; Friston et al., 2008), despite potentially failing to fully recover the spatial extent of the sources (Grova et al. 2006). The same was observed for inversion algorithms of the family of beamformers as the EBB used here, namely the dynamic imaging of coherent sources (DICS) and linearly constrained minimum variance (LCMV), exhibiting higher focal activation and lower spatial extent than those of LORETA (Halder et al. 2019). Interestingly, our post-hoc interaction analyses showed that by combining EBB with S2 or MSP with S3, the best performance in terms of d and p, respectively, is achieved, suggesting that by coupling inversion algorithms designed for providing focal solutions with information derived from fMRI data, an improved balance between specificity and sensitivity can be found. Furthermore, this type of comparison studies usually doesn’t reveal the “ground truth”, thus using simulated data can be useful to establish the appropriateness and optimality of the different methods. However, here we consider having a proxy for the ground truth accounting for the behavior of real data. fMRI priors would help improve EEG source imaging results if and where the fMRI activation and EEG sources are consistent. These are the anatomical ROIs where task-specific activations are expected but do not depend on the data being used for the reconstruction. Localizing in the brain the sources responsible for generating scalp EEG signals has been critical for determining the underpinnings of brain function in general, and those associated with multiple neurological disorders. Here, we compared for the first time EEG source reconstruction methods within a group of healthy subjects and MS patients.

Reconstructing Sources from EEG Data Collected Simultaneously with fMRI

The accurate reconstruction of EEG sources is not only related with the appropriateness of the inversion models used, but also with the overall quality of the EEG signal (Liu et al. 2018). EEG data simultaneously acquired with fMRI is known to suffer from severe artifact contamination (Abreu et al. 2018), but state-of-the-art pre-processing pipelines as the one used here can now bring data quality to sufficiently high levels. Despite the potentially inevitable loss in data quality relative to EEG collected outside the MR scanner, the feasibility of reconstructing sources of EEG data acquired simultaneously with fMRI has already been demonstrated (Groening et al. 2009; Vulliemoz et al. 2009, 2010a, b; Siniatchkin et al. 2010), particularly using EEG caps with a conventional spatial coverage (32 or 64 channels) as the one used in this study. Moreover, a direct relationship between EEG sources and fMRI networks has already been established first for data collected separately (Liu et al. 2017), and then validated on data collected simultaneously (Abreu et al. 2020b), supporting the feasibility of these procedures on this more challenging scenario. More importantly, analyzing EEG and fMRI data collected simultaneously is especially critical when studying spontaneous brain activity as the one associated with RSNs, for instance (Abreu et al. 2018).

Another important aspect is that fMRI and EEG record very different neuronal activity, so it is reasonable to think that fMRI-derived spatial priors will likely have an impact on source modelled neuronal activity (Goense et al. 2012). Also, fMRI is thought to reflect well high-frequency gamma activity but to be anticorrelated with alpha oscillations so one may think that using fMRI as prior will thus lead to a biased view where high-frequency activity is emphasized, and low-frequency activity dampened (Goense et al. 2012; Magri et al. 2012). Nevertheless, mapping high-frequency EEG oscillations with fMRI is not recommended due to contamination from residual artifacts, and in a worst-case scenario gamma activity could simply be muscle activity (Muthukumaraswamy 2013). Here, we addressed this issue by filtering out frequencies higher than 30 Hz from the EEG spectrum.

This further motivates the procedures performed here and suggest that deriving spatial priors from fMRI data separately acquired from EEG data may be suboptimal, which in turn could scale down their potential for guiding the reconstruction of EEG. Future studies would need to be conducted to confirm this observation.

Spatial Priors and Their Relationship with EEG (Sources)

In agreement with previous literature, in this work we found that including fMRI task activation maps and RSNs as additional CCs in the inversion models yielded significantly better EEG source reconstructions. This observation may be easily explained by the already known relationship between EEG and fMRI task activation and resting-state networks. In fact, source-reconstructed EEG data has already been used for mapping task-related fluctuations (Custo et al. 2014; Gonçalves et al. 2014), as well as for identifying RSNs (Liu et al. 2017, 2018; Abreu et al. 2020b), with recent studies showing a substantial overlap between these EEG maps and those typically obtained from fMRI data (Abreu et al. 2020b). Additionally, a relationship between fMRI RSNs and EEG has also been demonstrated in the sensor space, considering particularly the EEG rhythms extracted from the frequency domain (Goldman et al. 2002; Moosmann et al. 2003; Laufs et al. 2006; Scheeringa et al. 2008), which further supports the hypothesis that EEG carries in fact information that is also mapped with fMRI.

We then extended the exploration of fMRI spatial priors by also considering, for the first time, priors reflecting the fluctuations in the functional connectivity of task-related networks (dynamic functional connectivity, dFC). This was accomplished by estimating dFC fluctuations using phase coherence, followed by a dictionary learning step for finding the most recurrent dFC states, and a modularity analysis for identifying the network modules of the task-related dFC states. The rationale underlying our motivation for testing these spatial priors was based on recent literature showing that dFC fluctuations (Tagliazucchi et al. 2012; Chang et al. 2013; Preti et al. 2014; Korhonen et al. 2014; Tagliazucchi and Laufs 2015; Grooms et al. 2017; Omidvarnia et al. 2017), and dFC states in particular (Allen et al. 2018; Abreu et al. 2020a), have distinct EEG correlates, which could also be reflected on source reconstructed EEG data. Our results evidence this because by adding the task-related dFC state modules as spatial priors, the quality of the source reconstruction further increased for the task runs, in terms of the overlap with the ROIs and RSN templates. Noteworthy, other than task-related dFC states were not considered here, because otherwise all dFC states would have to be included given the lack of criteria for selecting a subset of them, which would be necessary to control for the potentially increasing complexity of the models.

Importantly, spatial priors of different natures have already been suggested (Lei et al. 2015). Knowing the structural connectome by analyzing diffusion MRI (dMRI) data may inform functional connectivity measures in the EEG source space in terms of the strength of the underlying structural connections, by weighting those measures accordingly (Knösche et al. 2013). Moreover, the fiber tracking between regions of interest allows to estimate the time lag between their functional connections, which is of particular interest when considering distantly located regions (Chu et al. 2015). Effective functional connectivity estimates obtained through Granger causality in the EEG source space can also be informed by connectivity priors also derived from Granger causality analyses of the fMRI data, despite its much lower temporal resolution compared to that of EEG (Roebroeck et al. 2005). Dynamic causal modeling (DCM) also estimates effective functional connectivity by incorporating information at the meso-scale (described by neural models whose parameters are typically defined based on animal studies) and the macro-scale (Friston et al. 2019). The latter has parameters reflecting forward, backward and lateral connections between sources, which can be defined from dMRI and/or structural MRI data. Similarly to Granger causality, DCM can also be applied to fMRI data, and the results used as connectivity priors for EEG source reconstruction (Lei et al. 2015). Naturally, these connectivity priors may be more crucial for studying EEG functional connectivity in the source space, which was not the case in the present study.

An alternative to the PEB framework for incorporating priors is the use of penalty functions (Lei et al. 2015). This constraint the inverse solutions using different types of norms and weight matrices that indirectly reflect a given prior, from which the MN and LORETA algorithms used here can be defined. Penalty functions have the advantage of easily balancing between sparse and smooth solutions by simply adjusting the norm accordingly, or combining multiple terms with different norms for intermediate solutions (Valdés-Sosa et al. 2009). However, the ability of explicitly incorporating spatial priors as covariance components in the inversion models designed under a PEB framework render it more interpretable, and therefore potentially more suitable for testing different types of fMRI spatial priors as it was performed here (López et al. 2014). Also, connectivity analysis in source space, which is a very promising approach in the context of MS, benefits from leveraging the advantages of both modalities. By leveraging simultaneous EEG-fMRI recordings, the incorporation of fMRI priors reveals its value by allowing to spatially and accurately locate brain activity captured by EEG, and to extract connectivity measures with a high temporal resolution. Notably, as MS is a disconnection disease in which the dynamics of the brain are altered, the information provided by the dFC state priors makes the reconstruction more reliable, by reflecting true neural mechanisms. On the other hand, the method might be useful from the point of view of the identification of unexpected sources, if these emerge from smaller signal generators that maybe associated with relatively weaker BOLD signal. Furthermore, it is also possible that unexpected sources arise from the solution of highly temporally resolved sources. This might be the case, e.g., if two spatially close neuronal populations would activate differentially in time in response to different task conditions, within the time scale of the EEG but not fMRI. However, if priors are not used, unexpected sources might be associated with a less reliable reconstruction, since EEG alone has poor spatial resolution to map brain task-activated regions. This yields a review of the applied methodology and a validation of the sources, particularly in cases in which the sources are unexpected regarding the task at hand. Therefore, this study also systematically compares different inversion algorithms without fMRI-priors while relying on the anatomical ROIs (where task-specific activations are expected) as a proxy for the ground truth. Thus, the investigation of the most appropriate or accurate inversion algorithm, even in the absence of fMRI priors in other studies, might be beneficial for the EEG/MEG connectivity community.

EEG Source Imaging in MS

The rationale for including patients with Multiple sclerosis (MS) in this study is that MS is a disconnection disease that is due to structural damage but also functional connectivity alterations. By leveraging the high temporal resolution from EEG and spatial resolution from fMRI, we can derive robust temporally- and spatially-resolved connectivity measures that might inform us better about the functional alterations. However, the reliability of the connectivity metrics depends on the reliability of the source reconstruction, in which fMRI-derived information can help improve source imaging solutions. Particularly, performing a task that elicits a specific known set of connected brain regions can highlight specific connectivity changes and thus help to better understand the pathophysiology of the disease. Moreover, more reliable, and informative connectivity measures, especially in longitudinal studies, might be the key to develop a tool for reliable disease progression assessment, which might improve the management of this condition and have a great impact in clinical needs. To reach this goal it is crucial to first test the kind of methodologies here presented, namely, to investigate the effect of the disease on the choice of method and set of priors.

The ANOVA test results revealed that the main effect of the group was not significant, which motivated us to proceed with the analyses with one group including all participants. Our rationale is that the quality of the sources’ reconstruction was irrespective of the presence of disease. We acknowledge the limited sample size; thus, these results should be seen as suggestive regarding the recommendations to future studies. As to our knowledge, there are no EEG-fMRI studies with focus on MS to identify connectivity biomarkers, so this work is a first and crucial step towards this goal. As MS can cause alterations to brain activity, which could influence the results, validation of these results is needed in future studies, namely in other healthy/patient cohorts alone with more data and considering other task designs. Nevertheless, this study represents a first step towards a more standard procedure for fMRI-informed EEG analyses in this context.

Conclusions

In this study, we systematically compared the quality of the source reconstruction of EEG data performed using different combinations of four inversion algorithms and three sets of covariance components incorporating different types of spatial priors derived from concurrently acquired fMRI data. We found that according to the quality metrics reflecting the presence of neuronal activity, combining the EBB or MSP algorithms with CC sets including fMRI task activation maps and RSNs yields the overall best source reconstruction, and that by further including dFC state modules as spatial priors, the quality of EEG sources from the task runs is optimal. We show that incorporating fMRI spatial priors in general, and for the first time dFC state modules in particular, plays a positive role in improving the reconstruction of EEG sources (and consequently any subsequent analyses). By providing a clear recommendation on the best approach for tackling the challenging inverse problem supported by our comprehensive analyses, we believe that future studies, particularly using real EEG data, may then be more informatively guided on this intricate research field.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 52 KB)