Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

EEG/MEG source imaging using fMRI informed time‐variant constraints

EEG/MEG source imaging using fMRI informed time‐variant constraints INTRODUCTIONElectroencephalography (EEG) and magnetoencephalography (MEG) can detect neural activities with high temporal resolution in milliseconds (Baillet, ). However, these two neuroimaging methods have a coarser spatial resolution due to the underdetermined inverse problems (Hämäläinen, Hari, Ilmoniemi, Knuutila, & Lounasmaa, ; Nunez and Srinivasan, ). Functional magnetic resonance imaging (fMRI), by measuring changes in blood oxygen level dependent (BOLD) signals, is able to provide a high spatial resolution on a scale of millimeters, but its temporal resolution is lower than that of MEG and EEG because of the slow hemodynamic responses (Kwong et al., ; Ogawa et al., ). Previous studies of neurovascular coupling suggested that BOLD response is linearly correlated with the power of local field potential (Logothetis, Pauls, Augath, Trinath, & Oeltermann, ; Logothetis and Wandell, ); thus, it is possible that a high spatial and temporal resolution mapping of brain activities can be generated by combining the information from EEG/MEG and fMRI (Dale and Halgren, ; Henson, Flandin, Friston, & Mattout, ; Liu, Belliveau, & Dale, ). Simulation and experimental evidence has already demonstrated the feasibility of combining EEG/MEG and fMRI data to obtain the mapping of brain activities with high spatiotemporal resolution (Dale et al., ; Henson et al., ; Liu, Dale, & Belliveau, ).However, fMRI constrained EEG or MEG source imaging relies on two assumptions. First, neural activities that can be detected by MEG or EEG are available in fMRI activation regions. Second, the neuronal activities are able to trigger vascular responses. As the dynamics of the neural and vascular activities differ substantially, and their exact relationship is yet to be fully characterized, situations that violate these two assumptions can be easily conceived (Ahlfors & Simpson, ). For instance, fMRI activation regions may remain electrically silent at the scalp due to the close‐field geometry (Huster, Debener, Eichele, & Herrmann, ). This example violates the first assumption, and the fMRI activation regions can be defined as fMRI extra sources. Under other circumstances, neural activities can be detected by EEG or MEG but cannot cause BOLD signal because the increase in cerebral blood flow caused by the activated neurons is too small to detect. This violates the second assumption and may cause fMRI missing sources. Due to the existence of fMRI extra sources or fMRI missing sources, fMRI constrained EEG or MEG source imaging may lead to inaccurate results (Fujimaki, Hayakawa, Nielsen, Knösche, & Miyauchi, ). Moreover, fMRI missing sources have greater negative impacts on the accuracy of fMRI constrained EEG or MEG source imaging (Liu et al., ; Liu and He, ).The approach that is used most often for fMRI constrained EEG or MEG source imaging is the fMRI‐weighted minimum norm estimation (fMNE) algorithm (Liu et al., ). To solve the mismatch problems between the location of fMRI activation regions and neural activities, the weights in fMNE are mainly based on fMRI activation maps. The diagonal elements of the source covariance matrix are set to 1 for the location in fMRI activation regions and 0.1 for other locations. The off‐diagonal elements are set to 0. They found that using this source covariance matrix provided the best tradeoff between sensitivities to crosstalk from fMRI visible and missing sources. In fMNE, those weights are assumed to be constant during the entire period of interest. However, neural activities change with time. Therefore, the constant weights used for all time points may cause excessive bias in the source waveforms of the EEG/MEG source estimation (Liu, Kecman, & He, ). Liu and He () proposed a time‐variant spatial constraints method which relies on the assumption that the BOLD signal is proportional to the time integral of local source power during the event related potential (ERP) generation.In this study, we propose a novel fMRI informed time‐variant constraint (FITC) method to address the weight problem mentioned above. Our proposed method does not rely on the assumption of the quantitative relation between the BOLD signal and neuronal activities. In contrast to fMNE, the weights in FITC are determined by both the neural electric weights and fMRI weights in a time‐variant manner to reduce the impact of fMRI extra sources. fMRI weights were then modified using crosstalk value ξ and normalized partial area under the curve value ρ to reduce the impact of fMRI missing sources. A location with lower ξ value represents higher estimation accuracy, and the location with a larger ρ value means that neurons are more active for a longer time. Thus, locations with lower ξ values and larger ρ values might represent the fMRI missing sources. “Depth weighting” of the lead field of FITC was applied to reduce the bias toward the superficial versus deep sources (Fuchs, Wagner, Kohler, & Wischmann, ; Lin et al., ).Computer simulations and experiments were performed to evaluate and compare the inverse solution based on the minimum norm estimation (MNE) (Hämäläinen and Ilmoniemi, ; Molins, Stufflebeam, Brown, & Hämäläinen, ), fMNE, FITC and depth weighted FITC (wFITC) method. In our simulation study, we randomly selected the dipole locations and simulated the corresponding fMRI activation region with/without introducing mismatches between fMRI activation and source locations to simulate the matched and mismatched conditions (fMRI missing source and extra source) at different signal‐to‐noise ratios (SNRs). Moreover, additional simulation factors (number of mismatched sources and temporally overlapping source waveforms) were also taken into account. Both the spatial and temporal accuracies were evaluated for conditions mentioned previously. In our visually evoked EEG/MEG and fMRI experiments, the source estimates responding to a visual stimulus was determined and compared among the following four algorithms: MNE, fMNE, FITC, and wFITC.METHODSEEG/MEG inverse problemThe forward model for EEG or MEG imaging can be expressed in a simple linear equation (Liu et al., ):1x(t)=Ls(t)+n(t)where x(t) is an Nx‐by‐1 vector (Nx is the total number of EEG electrodes or MEG sensors) of EEG or MEG recordings, s(t) is an Ns‐by‐1 vector (Ns is the total number of current sources) of unknown source strengths, L is an Nx‐by‐Ns gain matrix and n(t) is an Nx‐by‐1 noise vector.There are several ways to derive the linear inverse operator. The most frequently used are Tichonov regularization, Wiener filter, and Bayesian estimation (Hauk, ). All these derivations result in an equivalent linear inverse operator (G), which can be expressed as follows (Liu et al., ):2G= RLTLRLT+C−1where C is the noise covariance matrix and R is the source covariance matrix. Notice that if both the source covariance matrix R and noise covariance matrix C are set to the identity matrix, then the approach turns into the Moore–Penrose pseudoinverse (Hämäläinen and Ilmoniemi, 1994).Once the optimal linear inverse estimator G is calculated, the estimated source activity ŝ can be calculated using the simple expression:3ŝt=GxtfMRI‐weighted minimum norm estimationThe source covariance matrix or the weights R used in fMNE can be written as:4R=Rfwhere Rf is an Ns‐by‐Ns matrix. The diagonal elements of Rf are set to 1 for the locations in the fMRI activation regions, while other diagonal elements are set to 0.1. The off‐diagonal elements of Rf are set to 0 (Liu et al., ).fMRI informed time‐variant constraintsConstant weights for all time points may cause excess bias in the estimated source time courses in fMRI‐constrained EEG/MEG source estimation (Liu and He, ). To address this problem, the weights in our proposed FITC method are constructed in a time‐variant manner by multiplying the neural electric weights (i.e., Re) and fMRI weights (i.e., Rf):5R(t)=RfRe(t)Rf and Re are both diagonal matrices, Re can be constructed by taking its diagonal elements to be the MNE results, as follows:6 Re(t)=ŝ12(t)⋯0⋮ŝi2(t)⋮0⋯ŝNs2(t)where ŝi2 represents the estimate source power at location i. The weights in Equation highlight the fMRI activation regions with high source power. Then, the weights will reduce the effect of fMRI extra sources in fMRI‐constrained EEG/MEG source imaging.Because fMRI‐constrained EEG/MEG source imaging might result in inaccurate source estimates when there are fMRI missing sources, the fMRI weights Rf in Equation will be modified in the case of fMRI missing sources:7R(t)=Rξ, ρRe(t)Rξ, ρ is constructed based on two parameters: the normalized partial area under the curve matrix (i.e., ρi) and the crosstalk matrix (i.e., ξi) of MNE results (Liu et al., ).We defined a normalized partial area under the curve matrix ρi, as follows:8 ρi=∫0Tqi(t)∗ŝi2tdt∫0Tŝi2tdtwhere T is the entire time of interest, and qi(t) is an Ns‐by‐1 source spatial distribution mask vector with suprathreshold elements set to 1 and under‐threshold elements set to 0. The threshold of qi(t) is determined by t‐test between the baseline and task conditions for each location i. The statistical t‐value maps were corrected using Bonferroni correction at p < .05. The numerator of Equation is the sum of estimated activation power over the threshold at location i, while the denominator is the total power of the entire time of interest. A large ρi value indicates that the sources are more active for a longer time at location i.The crosstalk matrix ξi is defined as follows:9ξi=∑j≠iGiLj2GiLi2where Gi is the ith row of G, and Li is the ith column of L. The cross‐talk metric specifies the ratio of the contribution of the activity from all other locations j that is spread onto the location i to the contribution of activity at location i.ρave  and ξave are the averages of ρi and ξi over the intersection areas of the real fMRI activation regions and significant EEG/MEG response regions, respectively. Locations that are in significant EEG/MEG response regions but not in real fMRI activation regions and meet the conditions that ρi  is exceeding ρave and ξi is less than ξave will be defined as the extended activation regions. We set the diagonal elements of Rξ, ρ corresponding to real fMRI activation regions and the extended activation regions to 1; any other locations were set to 0.1. The off‐diagonal elements of Rξ, ρ were set to 0.The calculating process of the weights used in the FITC method is outlined in Figure . The weights address the problems related to both missing and extra fMRI sources.Schematic of the FITC weights calculation schemeComputer simulationsSimulations were applied to assess the performance of the FITC algorithm. Forward models were based on realistic‐geometry boundary element method (BEM) model reconstructed from the structural MRI images of three human subjects. The conductivities of the brain, skull and scalp were set to be 0.33, 0.0165, and 0.33 S/m, respectively (He et al., ; Lai et al., ; Oostendorp, Delbeke, & Stegeman, ). The source space consisted of 8196 vertices that were uniformly distributed on the cortical surface. Each vertex node is assumed to have one source oriented perpendicular to the surface. We assumed that the 64 electrodes were registered to the subject's scalp according to a 10–20 electrode system.The numbers of current dipoles and fMRI activation regions were altered to simulate five conditions that are typically encountered in human experiments: match, fMRI missing sources (one or two missing sources), and fMRI extra sources (one or two extra sources) conditions.Match conditionThree dipoles were placed on three selected locations of 8196 vertices, as shown in Figure a. The three dipole waveforms were simulated by quadratic function with the same peak amplitude and same duration of 100 ms. The peak latencies of the three dipole waveforms were set to be 50, 150, and 250 ms (as shown in the middle panel of Figure a). To meet the match condition, the fMRI‐activated areas were simulated to be centered at three selected points and all with a radius of 1 cm, as shown in the right panel of Figure a.Match condition (a): dipole locations (left), waveforms (middle), and fMRI responses (right); mismatch condition (b): fMRI missing source; mismatch condition (c): fMRI extra condition [Color figure can be viewed at wileyonlinelibrary.com]Mismatch condition: fMRI with one or two missing sourcesThe three dipole locations and waveforms were identical to those described in the match condition. The fMRI‐activated areas were simulated to be centered at one or two of three selected points, and all with a radius of 1 cm, leaving one or two dipoles to behave as fMRI missing sources, as shown in Figure b.Mismatch condition: fMRI with one or two extra sourcesThree locations were selected from the 8196 vertices. The three simulated fMRI‐activated areas were identical to those described in the match condition, as shown in the right panel of Figure c. Either one or two dipoles were placed at different selected locations, leaving one or two fMRI‐activated areas to behave as fMRI extra sources, as shown in Figure c. The dipole waveforms were simulated as a quadratic function with the same peak amplitude and duration of 100 ms. The peak latency is 250 ms under the two fMRI extra sources condition, whereas under the one fMRI extra source condition, the peak latencies were 50 and 250 ms for the “red” and “blue” sources, respectively.Under each condition, the simulated scalp EEG signals were produced by data of the BEM‐based forward model and Gaussian noise. The baseline period was defined as 200 ms prior to stimulus onset and was assumed to contain only Gaussian noise. EEG data sets were simulated with two different SNRs. The SNR was defined as the ratio between the maximum power of ideal EEG measurements and the standard deviation of the baseline period (Mantini, Franciotti, Romani, & Pizzella, ). The SNR of one data set was 10, which was within the typical SNR range of real EEG data. A low SNR (SNR = 3) was also applied to the simulation to assess the performance of the proposed method. Spatial dispersion (SD) was calculated to estimate the spatial extent for different algorithms (Hauk, Wakeman, & Henson, ; Molins et al., ).10SDi=∑jdijKij2∑jKij2where K is the “resolution matrix” and Kij = GiLj; dij is the distance between the locations j and i.The overlapping source waveforms setting was used to simulated EEG signals (as shown in the middle panel of Figure a). The amplitude and duration were identical to those described above. The peak latencies were set at 50, 80, and 110 ms. The SNR of the simulated EEG signal was 10, which was within the typical SNR range of real EEG data.Monte Carlo simulationThe above simulations were repeated 100 times for each condition and each head model to estimate the performance statistics of the four different algorithms (MNE, fMNE, FITC, and wFITC). For each simulation, the source locations were randomly selected. When comparing the four different algorithms, both the temporal and spatial accuracies of the source estimates of the “blue” source were evaluated. Temporal accuracy is defined as the correlation coefficient between the estimated source waveforms and preset source waveforms of the real dipole location. The localized error was defined as the distance between the real dipole location and the most active estimated source location at the time of peak latency. To test whether the performance of these four methods differed, one‐way Analysis of Variation (ANOVA) with Tukey's post‐hoc comparisons was applied using Matlab software (The MathWorks, Natick, USA).Visual stimulation experimentsTo validate the applicability of FITC, we collected human EEG, MEG, and fMRI data during a checkerboard visual stimulus experiment.Three healthy subjects were recruited to participate in the experiment. Written informed consent was obtained before scanning, and the protocol of this study was approved by the institutional review board of Peking University. A rectangular checkerboard was displayed in the lower left quadrant of the visual field to serve as the visual stimulation. In the EEG and MEG experiments, the visual stimuli were displayed for 0.5 s with an inter stimulus interval (uniform black field) of one second in each trial. There were 13 blocks in the experiment, and each block included 10 trials followed by five seconds of a uniform black field. In all conditions, a central fixation was present. Subjects were instructed to minimize eye blinks during the visual stimuli condition. In the fMRI experiment, there were fifteen 20 s blocks of visual stimuli with ten 20 s resting blocks in between.The MEG data were acquired using a 306‐channel Neuromag MEG system (Elekta, Helsinki, Finland) at the Center for MRI Research at Peking University, and the EEG data were simultaneously acquired with a 64‐channel, MEG‐compatible EEG system with a sampling rate of 1000 Hz. Both EEG and MEG data were filtered using a 0.5–40 Hz band‐pass filter and then segmented into epochs from 200 ms prior to stimulus onset, which lasted for 500 ms. Trials that contained ocular artifacts were rejected. After these preprocessing steps, all trials were averaged to obtain the visual evoked field (VEF) and potential (VEP).The structural and functional MRI data were collected on a Siemens Magnetom Prisma 3T MR system (Siemens Medical Solutions, Erlangen, Germany) at the Center for MRI Research at Peking University. The whole‐head structural MR images were acquired using the MPRAGE sequence (TR = 2,530 ms, TE = 2.98 ms, matrix size 448 × 512 × 192 and voxel size 0.5 × 0.5 × 1 mm3). The T2*‐weighted fMRI data were acquired using the echo planar imaging (EPI) sequence (TR = 2,000 ms, TE = 30 ms, matrix size 64 × 64 × 33 and voxel size 3.5 × 3.5 × 4.2 mm3). The structural MRI and fMRI data were analyzed using the SPM12 package (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/). The fMRI data underwent standard preprocessing steps, including slice timing correction, head motion correction (Friston, Williams, Howard, Frackowiak, & Turner, ), and spatially smoothing by convolving with an 8‐mm FWHM Gaussian kernel. The fMRI activation maps were compared via Student's t‐test (p ≤ .05, FWE‐corrected) comparing the fixation to checkerboard visual stimulation conditions (Friston, Jezzard, & Turner, ).Four algorithms (MNE, fMNE, FITC, and wFITC) were used to estimate the source distribution. The fMNE, FITC, and wFITC algorithms used both VEP/VEF and fMRI data, while the MNE algorithm used the VEP/VEF data alone.RESULTSSimulation studiesMatch conditionFigure shows the source estimates reconstructed from MNE, fMNE, FITC and wFITC algorithms at the three peak latencies (50, 150, and 250 ms) under the match condition. The MNE algorithm tended to produce wide spread areas of current density, which is common for EEG or MEG source reconstruction (first row of Figure ). The fMNE approach imposed false positive regions. FITC and wFITC provided a better estimate of the spatial extent of the activations than MNE and fMNE (third and fourth rows of Figure ).Source estimates under the match condition for one dataset. Source estimates at three peak latencies (50, 150, and 250 ms) using four different algorithms (MNE, fMNE, FITC, and wFITC). All images were scaled to their individual maximum. The value at the top of the bars indicates the relative maximum amplitude relative to the left column [Color figure can be viewed at wileyonlinelibrary.com]Mismatch condition: fMRI with one or two missing sourcesThe left panel of Figure shows the source estimates reconstructed from MNE, fMNE, FITC and wFITC at the peak latencies (50, 150, and 250 ms) under the one fMRI missing sources condition. Both FITC and wFITC provided better spatial extent of fMRI visible sources than MNE at 50 and 150 ms (first row of Figure ). However, the fMRI missing source (i.e., the “blue” source) was considerably underestimated using the fMNE algorithm (second row of Figure ), and the most active source location was in the fMRI activated region that centered in the “red” source. FITC and wFITC improved the estimation accuracy compared to that of MNE and fMNE under one fMRI missing source condition (third and fourth rows of Figure ).Source estimates under the fMRI missing sources condition for one dataset. Left panel: one fMRI missing source condition; Right panel: two fMRI missing sources condition; Source estimates at three peak latencies (50, 150, and 250 ms) using four different algorithms (MNE, fMNE, FITC, and wFITC). All images were scaled to their individual maximum. The value at the top of the bars indicates the relative maximum amplitude relative to the left column [Color figure can be viewed at wileyonlinelibrary.com]When there were two fMRI missing sources, fMNE, FITC, and wFITC provided better estimated spatial extent of the fMRI visible source compared with MNE at 50 ms. However, the fMRI missing sources at 150 and 250 ms were underestimated by using fMNE as shown in the right panel of Figure .Mismatch condition: fMRI with one or two extra sourcesThe left panel of Figure shows the source estimates reconstructed from MNE, fMNE, FITC and wFITC algorithms at the peak latencies (50, 150, and 250 ms) under the fMRI extra sources condition. All methods could localize the sources, as shown in Figure . fMNE displayed small false positive regions at fMRI extra source locations. When the fMRI active regions covered the dipoles, both the FITC and wFITC algorithms yielded more narrow spread source images around pre‐set source locations than the MNE solution in Figure .Source estimates under the fMRI extra source condition for one dataset. Left panel: one fMRI extra source condition; right panel: two fMRI extra sources condition. Source estimates at three peak latencies (50, 150, and 250 ms) using four different algorithms (MNE, fMNE, FITC, and wFITC). All images were scaled to their individual maximum. The value at the top of the bars indicates the relative maximum amplitude relative to the left column [Color figure can be viewed at wileyonlinelibrary.com]When there were two fMRI extra sources, fMNE displayed false‐positive regions at fMRI extra source locations. In contrast, FITC and wFITC suppressed the weights of the fMRI extra source regions and the results were more accurate, as shown in the right panel of Figure .Monte Carlo simulationWhen the SNR of the simulated EEG signals is 10, as shown in the first row of Figure , fMNE, FITC and wFITC displayed an increase in temporal correlation relative to the EEG‐alone MNE solution. fMNE, FITC, and wFITC displayed a significant decrease in localization error than MNE under the match condition.Mean value and standard error of localization error and temporal accuracy at two SNRs for the following four methods: MNE, fMNE, FITC, and wFITC. Left panels: localization error; right panels: temporal accuracy. One‐way ANOVA with Tukey's post‐hoc comparisons (*p < .05, **p < .01, ***p < .001) [Color figure can be viewed at wileyonlinelibrary.com]When there was one fMRI missing source as show in Figure , fMNE suffered the lowest temporal correlation and spatial accuracy among these four algorithms. MNE, FITC, and wFITC displayed a significant increase in temporal correlation and spatial accuracy compared with fMNE.When there were two fMRI missing sources, the localization error and the temporal accuracy displayed similar patterns with the one fMRI missing source condition. MNE, FITC, and wFITC displayed a significant increase in temporal correlation and spatial accuracy compared with fMNE.When there was one fMRI extra source, FITC and wFITC achieved higher temporal correlation than MNE and fMNE. fMNE and wFITC displayed a significant increase in spatial accuracy relative to MNE.When there were two fMRI extra sources, the localization error and the temporal accuracy of the four methods displayed a similar pattern as the one fMRI extra source condition.The temporal and spatial accuracies of lower SNR (SNR = 3) was displayed in the second row of Figure . All four method displayed lower temporal and spatial accuracies than under higher SNR (SNR = 10). FITC and wFITC also displayed a significant increase in temporal and spatial accuracies relative to MNE and fMNE under the fMRI missing sources condition. Under the match condition, FITC and wFITC displayed a significant increase in temporal and spatial accuracies relative to MNE. fMNE, FITC and wFITC achieved higher temporal and spatial accuracies compared with MNE under the fMRI extra sources condition.Under the overlapping source waveforms configuration, all the algorithms displayed lower temporal correlation and higher localization error relative to the nonoverlapping source waveform at the same SNR as show in Figure . Under the match and fMRI extra source conditions, FITC, wFITC, and fMNE achieved higher accuracy than MNE. Under the fMRI missing source conditions, FITC and wFITC displayed significant increases in temporal and spatial accuracies compared with fMNE and MNE.Mean value and standard error of localization error and temporal accuracy under the overlapping source waveforms setting for the following three methods: MNE, fMNE, FITC, and wFITC. Left panels: localization error; right panels: temporal accuracy. One‐way ANOVA with Tukey's post‐hoc comparisons (*p < .05, **p < .01, ***p < .001) [Color figure can be viewed at wileyonlinelibrary.com]Visual stimulation experimentsThe fMRI activation map (p ≤ .05, FWE‐corrected) corresponding to visual stimulation is presented in Figure a. Brodmann area 17 (V1), area 18 (V2), and area 19 (V3, V4, V5) are shown in Figure b. The first through fourth rows of Figure c show the source estimates at 92, 115, and 175 ms after stimulus onset using the four source reconstruction algorithms (MNE, fMNE, FITC, and wFITC) on inflated cortical surfaces, respectively. From the source estimates obtained from MNE, the activations gradually extended from the low‐tier to high‐tier visual areas (first row in Figure c). Using our proposed FITC and wFITC algorithms to reconstruct EEG source images, dynamic neural activities starting from V1 to V2, V3, V4, and V5 could be observed with enhanced spatial resolution (third and fourth rows in Figure c). In contrast, V5 is not present in the results obtained using the fMNE at the latency of 175 ms, indicating that the vascular activities in this area is too weak for fMRI measurement (second row in Figure c).fMRI activation map (a); Brodmann area (b); source estimates of EEG data (c) at three latencies (92, 115, and 175 ms after the visual onset) using MNE (first row), fMNE (second row), FITC (third row), and wFITC (fourth row) [Color figure can be viewed at wileyonlinelibrary.com]The fMRI activation map is present in Figure a. The first through fourth rows of Figure c show the reconstructed MEG source data at 109, 135, and 175 ms after stimulus onset using the four source reconstruction algorithms (MNE, fMNE, FITC, and wFITC) on inflated cortical surfaces, respectively. Activity occurred at 109 ms in V1 using MNE, which was followed by strong, widespread activation in the dorsal and ventral pathways at 135 and 175 ms (first row in Figure c). As shown in the second row in Figure c, fMNE had improved spatial resolution compared with that of the MNE and showed stronger weighting toward the fMRI. Brodmann area 19 is present in MNE, but not in fMNE at 175 ms. FITC and wFITC removed several diffuse activation areas in the MNE results and detected activation in the visual area (third and fourth rows in Figure c). The source estimates of the other two subjects were presented in Supporting Information, Figures S3–S6.fMRI activation map (a); Brodmann area (b); source estimates of MEG data (c) at three latencies (109, 135, and 175 ms after the visual onset) using MNE (first row), fMNE (second row), FITC (third row), and wFITC (fourth row) [Color figure can be viewed at wileyonlinelibrary.com]DISCUSSIONIn this study, we proposed an fMRI informed time‐variant constrained (FITC) algorithm to solve the mismatch problem between the locations of neural activities and fMRI activation regions. FITC weights were constructed in a time‐variant manner and determined by both EEG or MEG data and fMRI data. Monte Carlo simulation and visually evoked EEG, MEG and fMRI experiments were implemented to validate FITC. The simulation results demonstrate that the proposed FITC algorithm can handle fMRI missing and extra sources at different SNRs. The depth weighted function was also applied to FITC. The simulation results demonstrate that the temporal and spatial accuracies of FITC and wFITC were higher than MNE under both the match and mismatched (fMRI extra or missing sources) conditions. Moreover, the spatial and temporal accuracies of FITC were significantly higher than those of fMNE. The overlapping source waveforms settings were also implemented in the simulation. The results demonstrate that FITC and wFITC achieved higher spatial and temporal accuracies under all conditions compared to MNE. Moreover, FITC and wFITC achieved higher accuracies than fMNE under the fMRI missing sources condition.In the visually evoked experiments, the estimated source distribution produced by FITC had a narrow spread, and we observed dynamic activity from lower‐tier to high‐tier visual areas. This observed dynamic activity generally agrees with previous studies (Bonmassar, Schwartz, & Liu, ; Cottereau, Ales, & Norcia, ; Vanni et al., ). There are differences between the activation pattern of the source estimates of EEG and MEG data since EEG and MEG have different sensitivities to different configured brain activities. EEG can detect both radial and tangential components of the electromagnetic field, while MEG shows higher sensitivity to the tangential component though it is insensitive to the radial component. The experimental results also showed that FITC effectively estimated the source distribution of the MEG and EEG data and provides high spatial resolution of brain activities.The reason that the FITC weights were constructed in a time‐variant manner is that neural activities change with time in the total time period of interest. Neural activities at a certain time point may only involve a subset of the fMRI activation regions, so that other activation regions appear as fMRI extra sources (Liu et al., ; Ou et al., ). To highlight the fMRI activation regions with neural activities and to weaken the fMRI extra source, the FITC weights were constructed by multiplying the estimated source power and the fMRI weights. Therefore, FITC weights emphasize the fMRI activation regions with high source estimates. Moreover, compared to using only fMRI weights, applying FITC weights also weaken the negative impact of fMRI missing sources. To further weaken the impact of the fMRI missing sources, the fMRI weights were modified by the following two parameters: the cross‐talk effect and the normalized partial area under curve. A location with lower crosstalk value represents higher estimation accuracy, and the location with a larger ρ value means that neurons are more active for a longer time. Locations with lower crosstalk values and larger ρ values might represent the fMRI missing sources. Then, the weights of these locations are set to 1 to reduce the effect of fMRI missing sources.The simulation results reported above are based on three dipoles. We also investigated the source estimates of using five dipoles. The spatial and temporal accuracies is a slightly lower using the five dipole setting than using the three dipoles setting (Supporting Information, Figure S1). It also should be noted that the proposed FITC method depends on the initial spatial distribution of the MEG or EEG data. Standard MNE is a reasonable choice for initialization, since it is unbiased to fMRI (Fuchs et al., ; Lin et al., ). Point sources were also implemented in this study for the convenience of comparison. In future works, we plan to simulate EEG or MEG data using distributed sources.CONCLUSIONSThe proposed FITC method estimates the weights using both fMRI and EEG or MEG data in a time‐variant manner. Both the simulation and experimental results demonstrated that FITC led to more accurate source estimates than fMNE and MNE for various conditions, including the match, fMRI missing sources and fMRI extra sources conditions.ACKNOWLEDGMENTSThis work was supported by the National Key Research and Development Program of China (2017YFC0108900); National Natural Science Foundation of China grants (81430037, 81727808, 81790650, 81790651, 31530031, 81471376); Beijing Municipal Science & Technology Commission (Z171100000117012, Z161100000216152); Shenzhen Peacock Plan (KQTD2015033016104926); and Guangdong Pearl River Talents Plan Innovative and Entrepreneurial Team grant (2016ZT06S220). The authors thank National Center for Protein Sciences at Peking University in Beijing, China, for assistance with MRI data acquisition and data analyses.REFERENCESAhlfors, S. P., & Simpson, G. V. (2004). Geomatrixal interpretation of fMRI‐guided MEG/EEG inverse estimates. NeuroImage, 22, 323–332.Baillet, S. (2017). Magnetoencephalography for brain electrophysiology and imaging. Nature Neuroscience, 20, 327–339.Bonmassar, G., Schwartz, D. P., & Liu, A. K. (2001). Spatiotemporal brain imaging of visual‐evoked activity using interleaved EEG and fMRI recordings. NeuroImage, 13, 1035–1043.Cottereau, B. R., Ales, J. M., & Norcia, A. M. (2015). How to use fMRI functional localizers to improve EEG/MEG source estimation. Journal of Neuroscience Methods, 250, 64–73.Dale, A. M., & Halgren, E. (2001). Spatiotemporal mapping of brain activity by integration of multiple imaging modalities. Current Opinion in Neurobiology, 11, 202–208.Dale, A. M., Liu, A. K., Fischl, B. R., Buckner, R. L., Belliveau, J. W., Lewine, J. D., & Halgren, E. (2000). Dynamic statistical paramatrix mapping: Combining fMRI and MEG for high‐resolution imaging of cortical activity. Neuron, 26, 55–67.Friston, K. J., Jezzard, P., & Turner, R. (1994). Analysis of functional MRI time‐series. Human Brain Mapping, 1, 153–171.Friston, K. J., Williams, S., Howard, R., Frackowiak, R. S. J., & Turner, R. (1996). Movement‐related effects in fMRI time‐series. Magnetic Resonance in Medicine, 35, 346–355.Fujimaki, N., Hayakawa, T., Nielsen, M., Knösche, T. R., & Miyauchi, S. (2002). An fMRI‐constrained MEG source analysis with procedures for dividing and grouping activation. NeuroImage, 17, 324–343.Fuchs, M., Wagner, M., Kohler, T., & Wischmann, H. A. (1999). Linear and nonlinear current density reconstructions. Journal of Clinical Neurophysiology, 16, 267–295.Hämäläinen, M. S., & Ilmoniemi, R. J. (1984). Interpreting measured magnetic fields of the brain: estimates of current distributions. Helsinki University of Technology, Department of Technical Physics Report TKK‐F‐A559.Hämäläinen, M. S., Hari, R., Ilmoniemi, R. J., Knuutila, J., & Lounasmaa, O. V. (1993). Magnetoencephalography—theory, instrumentation, and applications to noninvasive studies of the working human brain. Reviews of Modern Physics, 65, 413–497.Hämäläinen, M. S., & Ilmoniemi, R. J. (1994). Interpreting magnetic fields of the brain: Minimum norm estimates. Medical &Amp; Biological Engineering &Amp; Computing, 32, 35–42.Hauk, O. (2004). Keep it simple: A case for using classical minimum norm estimation in the analysis of EEG and MEG data. NeuroImage, 21, 1612–1621.Hauk, O., Wakeman, D. G., & Henson, R. N. (2011). Comparison of noise‐normalized minimum norm estimates for MEG analysis using multiple resolution metrics. NeuroImage, 54, 1966–1974.He, B., Musha, T., Okamoto, Y., Homma, S., Nakajima, Y., & Sato, T. (1987). Electric dipole tracing in the brain by means of the boundary element method and its accuracy. IEEE Transactions on Bio‐Medical Engineering, 34, 406–414.Henson, R. N., Flandin, G., Friston, K. J., & Mattout, J. (2010). A parametric empirical Bayesian framework for fMRI‐constrained MEG/EEG source reconstruction. Human Brain Mapping, 31, 1512–1531.Henson, R. N., Goshen‐Gottstein, Y., Ganel, T., Otten, L. J., Quayle, A., & Rugg, M. D. (2003). Electrophysiological and haemodynamic correlates of face perception, recognition and priming. Cerebral Cortex (New York, N.Y. : 1991), 13, 793–805.Huster, R. J., Debener, S., Eichele, T., & Herrmann, C. S. (2012). Methods for simultaneous EEG‐fMRI: An introductory review. The Journal of Neuroscience, 32, 6053–6060.Kwong, K. K., Belliveau, J. W., Chesler, D. A., Goldberg, I. E., Weisskoff, R. M., Poncelet, B. P., … Turner, R. (1992). Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation. Proceedings of the National Academy of Sciences of the United States of America, 89, 5675–5679.Lai, Y., Van Drongelen, W., Ding, L., Hecox, K., Towle, V., Frim, D., & He, B. (2005). Estimation of in vivo human brain‐to‐skull conductivity ratio from simultaneous extra‐and intra‐cranial electrical potential recordings. Clinical Neurophysiology, 116, 456–465.Lin, F. H., Witzel, T., Ahlfors, S. P., Stufflebeam, S. M., Belliveau, J. W., & Hämäläinen, M. S. (2006). Assessing and improving the spatial accuracy in MEG source localization by depth‐weighted minimum‐norm estimates. NeuroImage, 31, 160–171.Liu, A. K., Belliveau, J. W., & Dale, A. M. (1998). Spatiotemporal imaging of human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations. Proceedings of the National Academy of Sciences of the United States of America, 95, 8945–8950.Liu, A. K., Dale, A. M., & Belliveau, J. W. (2002). Monte Carlo simulation studies of EEG and MEG localization accuracy. Human Brain Mapping, 16, 47–62.Liu, Z., Kecman, F., & He, B. (2006). Effects of fMRI–EEG mismatches in cortical current density estimation integrating fMRI and EEG: A simulation study. Clinical Neurophysiology, 117, 1610–1622.Liu, Z., & He, B. (2008). fMRI–EEG integrated cortical source imaging by use of time‐variant spatial constraints. NeuroImage, 39, 1198–1214.Logothetis, N. K., Pauls, J., Augath, M., Trinath, T., & Oeltermann, A. (2001). Neurophysiological investigation of the basis of the fMRI signal. Nature, 412, 150–157.Logothetis, N. K., & Wandell, B. A. (2004). Interpreting the BOLD signal. Annual Review of Physiology, 66, 735–769.Mantini, D., Franciotti, R., Romani, G. L., & Pizzella, V. (2008). Improving MEG source localizations: An automated method for complete artifact removal based on independent component analysis. NeuroImage, 40, 160–173.Molins, A., Stufflebeam, S. M., Brown, E. N., & Hämäläinen, M. S. (2008). Quantification of the benefit from integrating MEG and EEG data in minimum ℓ2‐norm estimation. NeuroImage, 42, 1069–1077.Nunez, P., & Srinivasan, R. (2006). Electrical fields of the brain: The neurophysics of EEG, 2nd ed. Oxford University Press. p. 106.Ogawa, S., Tank, D. W., Menon, R., Ellermann, J. M., Kim, S. G., Merkle, H., & Ugurbil, K. (1992). Intrinsic signal changes accompanying sensory stimulation: Functional brain mapping with magnetic resonance imaging. Proceedings of the National Academy of Sciences of the United States of America, 89, 5951–5955.Oostendorp, T. F., Delbeke, J., & Stegeman, D. F. (2000). The conductivity of the human skull: Results of in vivo and in vitro measurements. IEEE Transactions on Bio‐Medical Engineering, 47, 1487–1492.Ou, W., Nummenmaa, A., Ahveninen, J., Belliveau, J. W., Hämäläinen, M. S., & Golland, P. (2010). Multimodal functional imaging using fMRI‐informed regional EEG/MEG source estimation. NeuroImage, 52, 97–108.Vanni, S., Warnking, J., Dojat, M., Delon‐Martin, C., Bullier, J., & Segebarth, C. (2004). Sequence of pattern onset responses in the human visual areas: An fMRI constrained VEP source analysis. NeuroImage, 21, 801–817. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Human Brain Mapping Wiley

EEG/MEG source imaging using fMRI informed time‐variant constraints

Download PDF
12 pages

Loading next page...
 
/lp/wiley/eeg-meg-source-imaging-using-fmri-informed-time-variant-constraints-RLaTqttODa

References (36)

Publisher
Wiley
Copyright
© 2018 Wiley Periodicals, Inc.
ISSN
1065-9471
eISSN
1097-0193
DOI
10.1002/hbm.23945
Publisher site
See Article on Publisher Site

Abstract

INTRODUCTIONElectroencephalography (EEG) and magnetoencephalography (MEG) can detect neural activities with high temporal resolution in milliseconds (Baillet, ). However, these two neuroimaging methods have a coarser spatial resolution due to the underdetermined inverse problems (Hämäläinen, Hari, Ilmoniemi, Knuutila, & Lounasmaa, ; Nunez and Srinivasan, ). Functional magnetic resonance imaging (fMRI), by measuring changes in blood oxygen level dependent (BOLD) signals, is able to provide a high spatial resolution on a scale of millimeters, but its temporal resolution is lower than that of MEG and EEG because of the slow hemodynamic responses (Kwong et al., ; Ogawa et al., ). Previous studies of neurovascular coupling suggested that BOLD response is linearly correlated with the power of local field potential (Logothetis, Pauls, Augath, Trinath, & Oeltermann, ; Logothetis and Wandell, ); thus, it is possible that a high spatial and temporal resolution mapping of brain activities can be generated by combining the information from EEG/MEG and fMRI (Dale and Halgren, ; Henson, Flandin, Friston, & Mattout, ; Liu, Belliveau, & Dale, ). Simulation and experimental evidence has already demonstrated the feasibility of combining EEG/MEG and fMRI data to obtain the mapping of brain activities with high spatiotemporal resolution (Dale et al., ; Henson et al., ; Liu, Dale, & Belliveau, ).However, fMRI constrained EEG or MEG source imaging relies on two assumptions. First, neural activities that can be detected by MEG or EEG are available in fMRI activation regions. Second, the neuronal activities are able to trigger vascular responses. As the dynamics of the neural and vascular activities differ substantially, and their exact relationship is yet to be fully characterized, situations that violate these two assumptions can be easily conceived (Ahlfors & Simpson, ). For instance, fMRI activation regions may remain electrically silent at the scalp due to the close‐field geometry (Huster, Debener, Eichele, & Herrmann, ). This example violates the first assumption, and the fMRI activation regions can be defined as fMRI extra sources. Under other circumstances, neural activities can be detected by EEG or MEG but cannot cause BOLD signal because the increase in cerebral blood flow caused by the activated neurons is too small to detect. This violates the second assumption and may cause fMRI missing sources. Due to the existence of fMRI extra sources or fMRI missing sources, fMRI constrained EEG or MEG source imaging may lead to inaccurate results (Fujimaki, Hayakawa, Nielsen, Knösche, & Miyauchi, ). Moreover, fMRI missing sources have greater negative impacts on the accuracy of fMRI constrained EEG or MEG source imaging (Liu et al., ; Liu and He, ).The approach that is used most often for fMRI constrained EEG or MEG source imaging is the fMRI‐weighted minimum norm estimation (fMNE) algorithm (Liu et al., ). To solve the mismatch problems between the location of fMRI activation regions and neural activities, the weights in fMNE are mainly based on fMRI activation maps. The diagonal elements of the source covariance matrix are set to 1 for the location in fMRI activation regions and 0.1 for other locations. The off‐diagonal elements are set to 0. They found that using this source covariance matrix provided the best tradeoff between sensitivities to crosstalk from fMRI visible and missing sources. In fMNE, those weights are assumed to be constant during the entire period of interest. However, neural activities change with time. Therefore, the constant weights used for all time points may cause excessive bias in the source waveforms of the EEG/MEG source estimation (Liu, Kecman, & He, ). Liu and He () proposed a time‐variant spatial constraints method which relies on the assumption that the BOLD signal is proportional to the time integral of local source power during the event related potential (ERP) generation.In this study, we propose a novel fMRI informed time‐variant constraint (FITC) method to address the weight problem mentioned above. Our proposed method does not rely on the assumption of the quantitative relation between the BOLD signal and neuronal activities. In contrast to fMNE, the weights in FITC are determined by both the neural electric weights and fMRI weights in a time‐variant manner to reduce the impact of fMRI extra sources. fMRI weights were then modified using crosstalk value ξ and normalized partial area under the curve value ρ to reduce the impact of fMRI missing sources. A location with lower ξ value represents higher estimation accuracy, and the location with a larger ρ value means that neurons are more active for a longer time. Thus, locations with lower ξ values and larger ρ values might represent the fMRI missing sources. “Depth weighting” of the lead field of FITC was applied to reduce the bias toward the superficial versus deep sources (Fuchs, Wagner, Kohler, & Wischmann, ; Lin et al., ).Computer simulations and experiments were performed to evaluate and compare the inverse solution based on the minimum norm estimation (MNE) (Hämäläinen and Ilmoniemi, ; Molins, Stufflebeam, Brown, & Hämäläinen, ), fMNE, FITC and depth weighted FITC (wFITC) method. In our simulation study, we randomly selected the dipole locations and simulated the corresponding fMRI activation region with/without introducing mismatches between fMRI activation and source locations to simulate the matched and mismatched conditions (fMRI missing source and extra source) at different signal‐to‐noise ratios (SNRs). Moreover, additional simulation factors (number of mismatched sources and temporally overlapping source waveforms) were also taken into account. Both the spatial and temporal accuracies were evaluated for conditions mentioned previously. In our visually evoked EEG/MEG and fMRI experiments, the source estimates responding to a visual stimulus was determined and compared among the following four algorithms: MNE, fMNE, FITC, and wFITC.METHODSEEG/MEG inverse problemThe forward model for EEG or MEG imaging can be expressed in a simple linear equation (Liu et al., ):1x(t)=Ls(t)+n(t)where x(t) is an Nx‐by‐1 vector (Nx is the total number of EEG electrodes or MEG sensors) of EEG or MEG recordings, s(t) is an Ns‐by‐1 vector (Ns is the total number of current sources) of unknown source strengths, L is an Nx‐by‐Ns gain matrix and n(t) is an Nx‐by‐1 noise vector.There are several ways to derive the linear inverse operator. The most frequently used are Tichonov regularization, Wiener filter, and Bayesian estimation (Hauk, ). All these derivations result in an equivalent linear inverse operator (G), which can be expressed as follows (Liu et al., ):2G= RLTLRLT+C−1where C is the noise covariance matrix and R is the source covariance matrix. Notice that if both the source covariance matrix R and noise covariance matrix C are set to the identity matrix, then the approach turns into the Moore–Penrose pseudoinverse (Hämäläinen and Ilmoniemi, 1994).Once the optimal linear inverse estimator G is calculated, the estimated source activity ŝ can be calculated using the simple expression:3ŝt=GxtfMRI‐weighted minimum norm estimationThe source covariance matrix or the weights R used in fMNE can be written as:4R=Rfwhere Rf is an Ns‐by‐Ns matrix. The diagonal elements of Rf are set to 1 for the locations in the fMRI activation regions, while other diagonal elements are set to 0.1. The off‐diagonal elements of Rf are set to 0 (Liu et al., ).fMRI informed time‐variant constraintsConstant weights for all time points may cause excess bias in the estimated source time courses in fMRI‐constrained EEG/MEG source estimation (Liu and He, ). To address this problem, the weights in our proposed FITC method are constructed in a time‐variant manner by multiplying the neural electric weights (i.e., Re) and fMRI weights (i.e., Rf):5R(t)=RfRe(t)Rf and Re are both diagonal matrices, Re can be constructed by taking its diagonal elements to be the MNE results, as follows:6 Re(t)=ŝ12(t)⋯0⋮ŝi2(t)⋮0⋯ŝNs2(t)where ŝi2 represents the estimate source power at location i. The weights in Equation highlight the fMRI activation regions with high source power. Then, the weights will reduce the effect of fMRI extra sources in fMRI‐constrained EEG/MEG source imaging.Because fMRI‐constrained EEG/MEG source imaging might result in inaccurate source estimates when there are fMRI missing sources, the fMRI weights Rf in Equation will be modified in the case of fMRI missing sources:7R(t)=Rξ, ρRe(t)Rξ, ρ is constructed based on two parameters: the normalized partial area under the curve matrix (i.e., ρi) and the crosstalk matrix (i.e., ξi) of MNE results (Liu et al., ).We defined a normalized partial area under the curve matrix ρi, as follows:8 ρi=∫0Tqi(t)∗ŝi2tdt∫0Tŝi2tdtwhere T is the entire time of interest, and qi(t) is an Ns‐by‐1 source spatial distribution mask vector with suprathreshold elements set to 1 and under‐threshold elements set to 0. The threshold of qi(t) is determined by t‐test between the baseline and task conditions for each location i. The statistical t‐value maps were corrected using Bonferroni correction at p < .05. The numerator of Equation is the sum of estimated activation power over the threshold at location i, while the denominator is the total power of the entire time of interest. A large ρi value indicates that the sources are more active for a longer time at location i.The crosstalk matrix ξi is defined as follows:9ξi=∑j≠iGiLj2GiLi2where Gi is the ith row of G, and Li is the ith column of L. The cross‐talk metric specifies the ratio of the contribution of the activity from all other locations j that is spread onto the location i to the contribution of activity at location i.ρave  and ξave are the averages of ρi and ξi over the intersection areas of the real fMRI activation regions and significant EEG/MEG response regions, respectively. Locations that are in significant EEG/MEG response regions but not in real fMRI activation regions and meet the conditions that ρi  is exceeding ρave and ξi is less than ξave will be defined as the extended activation regions. We set the diagonal elements of Rξ, ρ corresponding to real fMRI activation regions and the extended activation regions to 1; any other locations were set to 0.1. The off‐diagonal elements of Rξ, ρ were set to 0.The calculating process of the weights used in the FITC method is outlined in Figure . The weights address the problems related to both missing and extra fMRI sources.Schematic of the FITC weights calculation schemeComputer simulationsSimulations were applied to assess the performance of the FITC algorithm. Forward models were based on realistic‐geometry boundary element method (BEM) model reconstructed from the structural MRI images of three human subjects. The conductivities of the brain, skull and scalp were set to be 0.33, 0.0165, and 0.33 S/m, respectively (He et al., ; Lai et al., ; Oostendorp, Delbeke, & Stegeman, ). The source space consisted of 8196 vertices that were uniformly distributed on the cortical surface. Each vertex node is assumed to have one source oriented perpendicular to the surface. We assumed that the 64 electrodes were registered to the subject's scalp according to a 10–20 electrode system.The numbers of current dipoles and fMRI activation regions were altered to simulate five conditions that are typically encountered in human experiments: match, fMRI missing sources (one or two missing sources), and fMRI extra sources (one or two extra sources) conditions.Match conditionThree dipoles were placed on three selected locations of 8196 vertices, as shown in Figure a. The three dipole waveforms were simulated by quadratic function with the same peak amplitude and same duration of 100 ms. The peak latencies of the three dipole waveforms were set to be 50, 150, and 250 ms (as shown in the middle panel of Figure a). To meet the match condition, the fMRI‐activated areas were simulated to be centered at three selected points and all with a radius of 1 cm, as shown in the right panel of Figure a.Match condition (a): dipole locations (left), waveforms (middle), and fMRI responses (right); mismatch condition (b): fMRI missing source; mismatch condition (c): fMRI extra condition [Color figure can be viewed at wileyonlinelibrary.com]Mismatch condition: fMRI with one or two missing sourcesThe three dipole locations and waveforms were identical to those described in the match condition. The fMRI‐activated areas were simulated to be centered at one or two of three selected points, and all with a radius of 1 cm, leaving one or two dipoles to behave as fMRI missing sources, as shown in Figure b.Mismatch condition: fMRI with one or two extra sourcesThree locations were selected from the 8196 vertices. The three simulated fMRI‐activated areas were identical to those described in the match condition, as shown in the right panel of Figure c. Either one or two dipoles were placed at different selected locations, leaving one or two fMRI‐activated areas to behave as fMRI extra sources, as shown in Figure c. The dipole waveforms were simulated as a quadratic function with the same peak amplitude and duration of 100 ms. The peak latency is 250 ms under the two fMRI extra sources condition, whereas under the one fMRI extra source condition, the peak latencies were 50 and 250 ms for the “red” and “blue” sources, respectively.Under each condition, the simulated scalp EEG signals were produced by data of the BEM‐based forward model and Gaussian noise. The baseline period was defined as 200 ms prior to stimulus onset and was assumed to contain only Gaussian noise. EEG data sets were simulated with two different SNRs. The SNR was defined as the ratio between the maximum power of ideal EEG measurements and the standard deviation of the baseline period (Mantini, Franciotti, Romani, & Pizzella, ). The SNR of one data set was 10, which was within the typical SNR range of real EEG data. A low SNR (SNR = 3) was also applied to the simulation to assess the performance of the proposed method. Spatial dispersion (SD) was calculated to estimate the spatial extent for different algorithms (Hauk, Wakeman, & Henson, ; Molins et al., ).10SDi=∑jdijKij2∑jKij2where K is the “resolution matrix” and Kij = GiLj; dij is the distance between the locations j and i.The overlapping source waveforms setting was used to simulated EEG signals (as shown in the middle panel of Figure a). The amplitude and duration were identical to those described above. The peak latencies were set at 50, 80, and 110 ms. The SNR of the simulated EEG signal was 10, which was within the typical SNR range of real EEG data.Monte Carlo simulationThe above simulations were repeated 100 times for each condition and each head model to estimate the performance statistics of the four different algorithms (MNE, fMNE, FITC, and wFITC). For each simulation, the source locations were randomly selected. When comparing the four different algorithms, both the temporal and spatial accuracies of the source estimates of the “blue” source were evaluated. Temporal accuracy is defined as the correlation coefficient between the estimated source waveforms and preset source waveforms of the real dipole location. The localized error was defined as the distance between the real dipole location and the most active estimated source location at the time of peak latency. To test whether the performance of these four methods differed, one‐way Analysis of Variation (ANOVA) with Tukey's post‐hoc comparisons was applied using Matlab software (The MathWorks, Natick, USA).Visual stimulation experimentsTo validate the applicability of FITC, we collected human EEG, MEG, and fMRI data during a checkerboard visual stimulus experiment.Three healthy subjects were recruited to participate in the experiment. Written informed consent was obtained before scanning, and the protocol of this study was approved by the institutional review board of Peking University. A rectangular checkerboard was displayed in the lower left quadrant of the visual field to serve as the visual stimulation. In the EEG and MEG experiments, the visual stimuli were displayed for 0.5 s with an inter stimulus interval (uniform black field) of one second in each trial. There were 13 blocks in the experiment, and each block included 10 trials followed by five seconds of a uniform black field. In all conditions, a central fixation was present. Subjects were instructed to minimize eye blinks during the visual stimuli condition. In the fMRI experiment, there were fifteen 20 s blocks of visual stimuli with ten 20 s resting blocks in between.The MEG data were acquired using a 306‐channel Neuromag MEG system (Elekta, Helsinki, Finland) at the Center for MRI Research at Peking University, and the EEG data were simultaneously acquired with a 64‐channel, MEG‐compatible EEG system with a sampling rate of 1000 Hz. Both EEG and MEG data were filtered using a 0.5–40 Hz band‐pass filter and then segmented into epochs from 200 ms prior to stimulus onset, which lasted for 500 ms. Trials that contained ocular artifacts were rejected. After these preprocessing steps, all trials were averaged to obtain the visual evoked field (VEF) and potential (VEP).The structural and functional MRI data were collected on a Siemens Magnetom Prisma 3T MR system (Siemens Medical Solutions, Erlangen, Germany) at the Center for MRI Research at Peking University. The whole‐head structural MR images were acquired using the MPRAGE sequence (TR = 2,530 ms, TE = 2.98 ms, matrix size 448 × 512 × 192 and voxel size 0.5 × 0.5 × 1 mm3). The T2*‐weighted fMRI data were acquired using the echo planar imaging (EPI) sequence (TR = 2,000 ms, TE = 30 ms, matrix size 64 × 64 × 33 and voxel size 3.5 × 3.5 × 4.2 mm3). The structural MRI and fMRI data were analyzed using the SPM12 package (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/). The fMRI data underwent standard preprocessing steps, including slice timing correction, head motion correction (Friston, Williams, Howard, Frackowiak, & Turner, ), and spatially smoothing by convolving with an 8‐mm FWHM Gaussian kernel. The fMRI activation maps were compared via Student's t‐test (p ≤ .05, FWE‐corrected) comparing the fixation to checkerboard visual stimulation conditions (Friston, Jezzard, & Turner, ).Four algorithms (MNE, fMNE, FITC, and wFITC) were used to estimate the source distribution. The fMNE, FITC, and wFITC algorithms used both VEP/VEF and fMRI data, while the MNE algorithm used the VEP/VEF data alone.RESULTSSimulation studiesMatch conditionFigure shows the source estimates reconstructed from MNE, fMNE, FITC and wFITC algorithms at the three peak latencies (50, 150, and 250 ms) under the match condition. The MNE algorithm tended to produce wide spread areas of current density, which is common for EEG or MEG source reconstruction (first row of Figure ). The fMNE approach imposed false positive regions. FITC and wFITC provided a better estimate of the spatial extent of the activations than MNE and fMNE (third and fourth rows of Figure ).Source estimates under the match condition for one dataset. Source estimates at three peak latencies (50, 150, and 250 ms) using four different algorithms (MNE, fMNE, FITC, and wFITC). All images were scaled to their individual maximum. The value at the top of the bars indicates the relative maximum amplitude relative to the left column [Color figure can be viewed at wileyonlinelibrary.com]Mismatch condition: fMRI with one or two missing sourcesThe left panel of Figure shows the source estimates reconstructed from MNE, fMNE, FITC and wFITC at the peak latencies (50, 150, and 250 ms) under the one fMRI missing sources condition. Both FITC and wFITC provided better spatial extent of fMRI visible sources than MNE at 50 and 150 ms (first row of Figure ). However, the fMRI missing source (i.e., the “blue” source) was considerably underestimated using the fMNE algorithm (second row of Figure ), and the most active source location was in the fMRI activated region that centered in the “red” source. FITC and wFITC improved the estimation accuracy compared to that of MNE and fMNE under one fMRI missing source condition (third and fourth rows of Figure ).Source estimates under the fMRI missing sources condition for one dataset. Left panel: one fMRI missing source condition; Right panel: two fMRI missing sources condition; Source estimates at three peak latencies (50, 150, and 250 ms) using four different algorithms (MNE, fMNE, FITC, and wFITC). All images were scaled to their individual maximum. The value at the top of the bars indicates the relative maximum amplitude relative to the left column [Color figure can be viewed at wileyonlinelibrary.com]When there were two fMRI missing sources, fMNE, FITC, and wFITC provided better estimated spatial extent of the fMRI visible source compared with MNE at 50 ms. However, the fMRI missing sources at 150 and 250 ms were underestimated by using fMNE as shown in the right panel of Figure .Mismatch condition: fMRI with one or two extra sourcesThe left panel of Figure shows the source estimates reconstructed from MNE, fMNE, FITC and wFITC algorithms at the peak latencies (50, 150, and 250 ms) under the fMRI extra sources condition. All methods could localize the sources, as shown in Figure . fMNE displayed small false positive regions at fMRI extra source locations. When the fMRI active regions covered the dipoles, both the FITC and wFITC algorithms yielded more narrow spread source images around pre‐set source locations than the MNE solution in Figure .Source estimates under the fMRI extra source condition for one dataset. Left panel: one fMRI extra source condition; right panel: two fMRI extra sources condition. Source estimates at three peak latencies (50, 150, and 250 ms) using four different algorithms (MNE, fMNE, FITC, and wFITC). All images were scaled to their individual maximum. The value at the top of the bars indicates the relative maximum amplitude relative to the left column [Color figure can be viewed at wileyonlinelibrary.com]When there were two fMRI extra sources, fMNE displayed false‐positive regions at fMRI extra source locations. In contrast, FITC and wFITC suppressed the weights of the fMRI extra source regions and the results were more accurate, as shown in the right panel of Figure .Monte Carlo simulationWhen the SNR of the simulated EEG signals is 10, as shown in the first row of Figure , fMNE, FITC and wFITC displayed an increase in temporal correlation relative to the EEG‐alone MNE solution. fMNE, FITC, and wFITC displayed a significant decrease in localization error than MNE under the match condition.Mean value and standard error of localization error and temporal accuracy at two SNRs for the following four methods: MNE, fMNE, FITC, and wFITC. Left panels: localization error; right panels: temporal accuracy. One‐way ANOVA with Tukey's post‐hoc comparisons (*p < .05, **p < .01, ***p < .001) [Color figure can be viewed at wileyonlinelibrary.com]When there was one fMRI missing source as show in Figure , fMNE suffered the lowest temporal correlation and spatial accuracy among these four algorithms. MNE, FITC, and wFITC displayed a significant increase in temporal correlation and spatial accuracy compared with fMNE.When there were two fMRI missing sources, the localization error and the temporal accuracy displayed similar patterns with the one fMRI missing source condition. MNE, FITC, and wFITC displayed a significant increase in temporal correlation and spatial accuracy compared with fMNE.When there was one fMRI extra source, FITC and wFITC achieved higher temporal correlation than MNE and fMNE. fMNE and wFITC displayed a significant increase in spatial accuracy relative to MNE.When there were two fMRI extra sources, the localization error and the temporal accuracy of the four methods displayed a similar pattern as the one fMRI extra source condition.The temporal and spatial accuracies of lower SNR (SNR = 3) was displayed in the second row of Figure . All four method displayed lower temporal and spatial accuracies than under higher SNR (SNR = 10). FITC and wFITC also displayed a significant increase in temporal and spatial accuracies relative to MNE and fMNE under the fMRI missing sources condition. Under the match condition, FITC and wFITC displayed a significant increase in temporal and spatial accuracies relative to MNE. fMNE, FITC and wFITC achieved higher temporal and spatial accuracies compared with MNE under the fMRI extra sources condition.Under the overlapping source waveforms configuration, all the algorithms displayed lower temporal correlation and higher localization error relative to the nonoverlapping source waveform at the same SNR as show in Figure . Under the match and fMRI extra source conditions, FITC, wFITC, and fMNE achieved higher accuracy than MNE. Under the fMRI missing source conditions, FITC and wFITC displayed significant increases in temporal and spatial accuracies compared with fMNE and MNE.Mean value and standard error of localization error and temporal accuracy under the overlapping source waveforms setting for the following three methods: MNE, fMNE, FITC, and wFITC. Left panels: localization error; right panels: temporal accuracy. One‐way ANOVA with Tukey's post‐hoc comparisons (*p < .05, **p < .01, ***p < .001) [Color figure can be viewed at wileyonlinelibrary.com]Visual stimulation experimentsThe fMRI activation map (p ≤ .05, FWE‐corrected) corresponding to visual stimulation is presented in Figure a. Brodmann area 17 (V1), area 18 (V2), and area 19 (V3, V4, V5) are shown in Figure b. The first through fourth rows of Figure c show the source estimates at 92, 115, and 175 ms after stimulus onset using the four source reconstruction algorithms (MNE, fMNE, FITC, and wFITC) on inflated cortical surfaces, respectively. From the source estimates obtained from MNE, the activations gradually extended from the low‐tier to high‐tier visual areas (first row in Figure c). Using our proposed FITC and wFITC algorithms to reconstruct EEG source images, dynamic neural activities starting from V1 to V2, V3, V4, and V5 could be observed with enhanced spatial resolution (third and fourth rows in Figure c). In contrast, V5 is not present in the results obtained using the fMNE at the latency of 175 ms, indicating that the vascular activities in this area is too weak for fMRI measurement (second row in Figure c).fMRI activation map (a); Brodmann area (b); source estimates of EEG data (c) at three latencies (92, 115, and 175 ms after the visual onset) using MNE (first row), fMNE (second row), FITC (third row), and wFITC (fourth row) [Color figure can be viewed at wileyonlinelibrary.com]The fMRI activation map is present in Figure a. The first through fourth rows of Figure c show the reconstructed MEG source data at 109, 135, and 175 ms after stimulus onset using the four source reconstruction algorithms (MNE, fMNE, FITC, and wFITC) on inflated cortical surfaces, respectively. Activity occurred at 109 ms in V1 using MNE, which was followed by strong, widespread activation in the dorsal and ventral pathways at 135 and 175 ms (first row in Figure c). As shown in the second row in Figure c, fMNE had improved spatial resolution compared with that of the MNE and showed stronger weighting toward the fMRI. Brodmann area 19 is present in MNE, but not in fMNE at 175 ms. FITC and wFITC removed several diffuse activation areas in the MNE results and detected activation in the visual area (third and fourth rows in Figure c). The source estimates of the other two subjects were presented in Supporting Information, Figures S3–S6.fMRI activation map (a); Brodmann area (b); source estimates of MEG data (c) at three latencies (109, 135, and 175 ms after the visual onset) using MNE (first row), fMNE (second row), FITC (third row), and wFITC (fourth row) [Color figure can be viewed at wileyonlinelibrary.com]DISCUSSIONIn this study, we proposed an fMRI informed time‐variant constrained (FITC) algorithm to solve the mismatch problem between the locations of neural activities and fMRI activation regions. FITC weights were constructed in a time‐variant manner and determined by both EEG or MEG data and fMRI data. Monte Carlo simulation and visually evoked EEG, MEG and fMRI experiments were implemented to validate FITC. The simulation results demonstrate that the proposed FITC algorithm can handle fMRI missing and extra sources at different SNRs. The depth weighted function was also applied to FITC. The simulation results demonstrate that the temporal and spatial accuracies of FITC and wFITC were higher than MNE under both the match and mismatched (fMRI extra or missing sources) conditions. Moreover, the spatial and temporal accuracies of FITC were significantly higher than those of fMNE. The overlapping source waveforms settings were also implemented in the simulation. The results demonstrate that FITC and wFITC achieved higher spatial and temporal accuracies under all conditions compared to MNE. Moreover, FITC and wFITC achieved higher accuracies than fMNE under the fMRI missing sources condition.In the visually evoked experiments, the estimated source distribution produced by FITC had a narrow spread, and we observed dynamic activity from lower‐tier to high‐tier visual areas. This observed dynamic activity generally agrees with previous studies (Bonmassar, Schwartz, & Liu, ; Cottereau, Ales, & Norcia, ; Vanni et al., ). There are differences between the activation pattern of the source estimates of EEG and MEG data since EEG and MEG have different sensitivities to different configured brain activities. EEG can detect both radial and tangential components of the electromagnetic field, while MEG shows higher sensitivity to the tangential component though it is insensitive to the radial component. The experimental results also showed that FITC effectively estimated the source distribution of the MEG and EEG data and provides high spatial resolution of brain activities.The reason that the FITC weights were constructed in a time‐variant manner is that neural activities change with time in the total time period of interest. Neural activities at a certain time point may only involve a subset of the fMRI activation regions, so that other activation regions appear as fMRI extra sources (Liu et al., ; Ou et al., ). To highlight the fMRI activation regions with neural activities and to weaken the fMRI extra source, the FITC weights were constructed by multiplying the estimated source power and the fMRI weights. Therefore, FITC weights emphasize the fMRI activation regions with high source estimates. Moreover, compared to using only fMRI weights, applying FITC weights also weaken the negative impact of fMRI missing sources. To further weaken the impact of the fMRI missing sources, the fMRI weights were modified by the following two parameters: the cross‐talk effect and the normalized partial area under curve. A location with lower crosstalk value represents higher estimation accuracy, and the location with a larger ρ value means that neurons are more active for a longer time. Locations with lower crosstalk values and larger ρ values might represent the fMRI missing sources. Then, the weights of these locations are set to 1 to reduce the effect of fMRI missing sources.The simulation results reported above are based on three dipoles. We also investigated the source estimates of using five dipoles. The spatial and temporal accuracies is a slightly lower using the five dipole setting than using the three dipoles setting (Supporting Information, Figure S1). It also should be noted that the proposed FITC method depends on the initial spatial distribution of the MEG or EEG data. Standard MNE is a reasonable choice for initialization, since it is unbiased to fMRI (Fuchs et al., ; Lin et al., ). Point sources were also implemented in this study for the convenience of comparison. In future works, we plan to simulate EEG or MEG data using distributed sources.CONCLUSIONSThe proposed FITC method estimates the weights using both fMRI and EEG or MEG data in a time‐variant manner. Both the simulation and experimental results demonstrated that FITC led to more accurate source estimates than fMNE and MNE for various conditions, including the match, fMRI missing sources and fMRI extra sources conditions.ACKNOWLEDGMENTSThis work was supported by the National Key Research and Development Program of China (2017YFC0108900); National Natural Science Foundation of China grants (81430037, 81727808, 81790650, 81790651, 31530031, 81471376); Beijing Municipal Science & Technology Commission (Z171100000117012, Z161100000216152); Shenzhen Peacock Plan (KQTD2015033016104926); and Guangdong Pearl River Talents Plan Innovative and Entrepreneurial Team grant (2016ZT06S220). The authors thank National Center for Protein Sciences at Peking University in Beijing, China, for assistance with MRI data acquisition and data analyses.REFERENCESAhlfors, S. P., & Simpson, G. V. (2004). Geomatrixal interpretation of fMRI‐guided MEG/EEG inverse estimates. NeuroImage, 22, 323–332.Baillet, S. (2017). Magnetoencephalography for brain electrophysiology and imaging. Nature Neuroscience, 20, 327–339.Bonmassar, G., Schwartz, D. P., & Liu, A. K. (2001). Spatiotemporal brain imaging of visual‐evoked activity using interleaved EEG and fMRI recordings. NeuroImage, 13, 1035–1043.Cottereau, B. R., Ales, J. M., & Norcia, A. M. (2015). How to use fMRI functional localizers to improve EEG/MEG source estimation. Journal of Neuroscience Methods, 250, 64–73.Dale, A. M., & Halgren, E. (2001). Spatiotemporal mapping of brain activity by integration of multiple imaging modalities. Current Opinion in Neurobiology, 11, 202–208.Dale, A. M., Liu, A. K., Fischl, B. R., Buckner, R. L., Belliveau, J. W., Lewine, J. D., & Halgren, E. (2000). Dynamic statistical paramatrix mapping: Combining fMRI and MEG for high‐resolution imaging of cortical activity. Neuron, 26, 55–67.Friston, K. J., Jezzard, P., & Turner, R. (1994). Analysis of functional MRI time‐series. Human Brain Mapping, 1, 153–171.Friston, K. J., Williams, S., Howard, R., Frackowiak, R. S. J., & Turner, R. (1996). Movement‐related effects in fMRI time‐series. Magnetic Resonance in Medicine, 35, 346–355.Fujimaki, N., Hayakawa, T., Nielsen, M., Knösche, T. R., & Miyauchi, S. (2002). An fMRI‐constrained MEG source analysis with procedures for dividing and grouping activation. NeuroImage, 17, 324–343.Fuchs, M., Wagner, M., Kohler, T., & Wischmann, H. A. (1999). Linear and nonlinear current density reconstructions. Journal of Clinical Neurophysiology, 16, 267–295.Hämäläinen, M. S., & Ilmoniemi, R. J. (1984). Interpreting measured magnetic fields of the brain: estimates of current distributions. Helsinki University of Technology, Department of Technical Physics Report TKK‐F‐A559.Hämäläinen, M. S., Hari, R., Ilmoniemi, R. J., Knuutila, J., & Lounasmaa, O. V. (1993). Magnetoencephalography—theory, instrumentation, and applications to noninvasive studies of the working human brain. Reviews of Modern Physics, 65, 413–497.Hämäläinen, M. S., & Ilmoniemi, R. J. (1994). Interpreting magnetic fields of the brain: Minimum norm estimates. Medical &Amp; Biological Engineering &Amp; Computing, 32, 35–42.Hauk, O. (2004). Keep it simple: A case for using classical minimum norm estimation in the analysis of EEG and MEG data. NeuroImage, 21, 1612–1621.Hauk, O., Wakeman, D. G., & Henson, R. N. (2011). Comparison of noise‐normalized minimum norm estimates for MEG analysis using multiple resolution metrics. NeuroImage, 54, 1966–1974.He, B., Musha, T., Okamoto, Y., Homma, S., Nakajima, Y., & Sato, T. (1987). Electric dipole tracing in the brain by means of the boundary element method and its accuracy. IEEE Transactions on Bio‐Medical Engineering, 34, 406–414.Henson, R. N., Flandin, G., Friston, K. J., & Mattout, J. (2010). A parametric empirical Bayesian framework for fMRI‐constrained MEG/EEG source reconstruction. Human Brain Mapping, 31, 1512–1531.Henson, R. N., Goshen‐Gottstein, Y., Ganel, T., Otten, L. J., Quayle, A., & Rugg, M. D. (2003). Electrophysiological and haemodynamic correlates of face perception, recognition and priming. Cerebral Cortex (New York, N.Y. : 1991), 13, 793–805.Huster, R. J., Debener, S., Eichele, T., & Herrmann, C. S. (2012). Methods for simultaneous EEG‐fMRI: An introductory review. The Journal of Neuroscience, 32, 6053–6060.Kwong, K. K., Belliveau, J. W., Chesler, D. A., Goldberg, I. E., Weisskoff, R. M., Poncelet, B. P., … Turner, R. (1992). Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation. Proceedings of the National Academy of Sciences of the United States of America, 89, 5675–5679.Lai, Y., Van Drongelen, W., Ding, L., Hecox, K., Towle, V., Frim, D., & He, B. (2005). Estimation of in vivo human brain‐to‐skull conductivity ratio from simultaneous extra‐and intra‐cranial electrical potential recordings. Clinical Neurophysiology, 116, 456–465.Lin, F. H., Witzel, T., Ahlfors, S. P., Stufflebeam, S. M., Belliveau, J. W., & Hämäläinen, M. S. (2006). Assessing and improving the spatial accuracy in MEG source localization by depth‐weighted minimum‐norm estimates. NeuroImage, 31, 160–171.Liu, A. K., Belliveau, J. W., & Dale, A. M. (1998). Spatiotemporal imaging of human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations. Proceedings of the National Academy of Sciences of the United States of America, 95, 8945–8950.Liu, A. K., Dale, A. M., & Belliveau, J. W. (2002). Monte Carlo simulation studies of EEG and MEG localization accuracy. Human Brain Mapping, 16, 47–62.Liu, Z., Kecman, F., & He, B. (2006). Effects of fMRI–EEG mismatches in cortical current density estimation integrating fMRI and EEG: A simulation study. Clinical Neurophysiology, 117, 1610–1622.Liu, Z., & He, B. (2008). fMRI–EEG integrated cortical source imaging by use of time‐variant spatial constraints. NeuroImage, 39, 1198–1214.Logothetis, N. K., Pauls, J., Augath, M., Trinath, T., & Oeltermann, A. (2001). Neurophysiological investigation of the basis of the fMRI signal. Nature, 412, 150–157.Logothetis, N. K., & Wandell, B. A. (2004). Interpreting the BOLD signal. Annual Review of Physiology, 66, 735–769.Mantini, D., Franciotti, R., Romani, G. L., & Pizzella, V. (2008). Improving MEG source localizations: An automated method for complete artifact removal based on independent component analysis. NeuroImage, 40, 160–173.Molins, A., Stufflebeam, S. M., Brown, E. N., & Hämäläinen, M. S. (2008). Quantification of the benefit from integrating MEG and EEG data in minimum ℓ2‐norm estimation. NeuroImage, 42, 1069–1077.Nunez, P., & Srinivasan, R. (2006). Electrical fields of the brain: The neurophysics of EEG, 2nd ed. Oxford University Press. p. 106.Ogawa, S., Tank, D. W., Menon, R., Ellermann, J. M., Kim, S. G., Merkle, H., & Ugurbil, K. (1992). Intrinsic signal changes accompanying sensory stimulation: Functional brain mapping with magnetic resonance imaging. Proceedings of the National Academy of Sciences of the United States of America, 89, 5951–5955.Oostendorp, T. F., Delbeke, J., & Stegeman, D. F. (2000). The conductivity of the human skull: Results of in vivo and in vitro measurements. IEEE Transactions on Bio‐Medical Engineering, 47, 1487–1492.Ou, W., Nummenmaa, A., Ahveninen, J., Belliveau, J. W., Hämäläinen, M. S., & Golland, P. (2010). Multimodal functional imaging using fMRI‐informed regional EEG/MEG source estimation. NeuroImage, 52, 97–108.Vanni, S., Warnking, J., Dojat, M., Delon‐Martin, C., Bullier, J., & Segebarth, C. (2004). Sequence of pattern onset responses in the human visual areas: An fMRI constrained VEP source analysis. NeuroImage, 21, 801–817.

Journal

Human Brain Mapping – Wiley

Published: Jan 1, 2018

Keywords: ; ; ; ;

There are no references for this article.