Extension of Quantifiable Modification of sLORETA for Induced Oscillatory Changes in Magnetoencephalography

Quantitative Multiple-Source Estimation Based on sLORETA-qm

We defined the magnetic field measured by the m-th sensor at time t as b_m(t), and a set of measured data as a column vector:

b(t)=[b₁(t),b₂(t),...,b_NC(t)]^T

where N_C is the total number of sensors and T indicates the matrix transpose. We defined the lead field vector from the i-th voxel to all sensors as I_i of N_C × 1 vector, and the composite lead field matrix for all voxels as

L = [I₁, I₂, …, I_i ,… I_Ni] (N_C × N_V matrix)

where N_V is the total number of voxels set in the brain area.

Using sLORETA, the estimated source on the i-th voxel ${\stackrel{}{s}}_i\left(t\right)$ is expressed as

where G is the gram matrix, defined as G = LL^T [4].

Although sLORETA has almost no localization bias, it has point spread function [3], and thus false-positive sources can be estimated around the actual source localization. As a result, we regarded the voxel with a spatially or locally maximum among adjacent voxels as the most reliable information for source estimation. From these perspectives, after source reconstruction by sLORETA, a voxel with spatially or locally maximum at each time point was chosen as a feasible source and used for further analyses.

To quantify current density estimated by sLORETA, we assumed the actual single source $s_k\left(t\right)$ to be on the k-th voxel and the lead field from $s_k\left(t\right)$ to sensors to be equal to I_k. Under these conditions, sensor signals from $s_k\left(t\right)$ are expressed as

Substituting Eq. (2) into Eq. (1), we obtain:

With regard to the voxel with spatially or locally maximum, the i-th voxel is considered to be close to the k-th voxel. Hence, Eq. (3) is expressed as

As G is constant after definition of brain area and l_k decreases in proportion to square of the distance between the source and sensors, Eq. (4) indicates that the deeper the location of the estimated sources, the smaller the current density that can be estimated as compared with the actual current density. We therefore quantified current density using the following equation to correct such errors:

where ${\stackrel{}{s}}_{kq}\left(t\right)$ is a quantitative value of ${\stackrel{}{s}}_k\left(t\right)$.

In cases involving multiple sources, Eq. (2) is described as

where N_S is the number of actual sources, $s_{k_j}\left(t\right)$ is the actual source on k_j voxel and $l_{k_j}$ is the lead field vector from $s_{k_j}\left(t\right)$ to sensors.

In the same manner as the derivation of Eq. (4), estimated sources are expressed as

where ${\stackrel{}{s}}_{k_1}\left(t\right),{\stackrel{}{s}}_{k_2}\left(t\right),...,{\stackrel{}{s}}_{k_{N_s}}\left(t\right)$ are estimated sources at $k_1,k_2,...,k_{N_s}$-th voxels. Finally, we can obtain quantified current density using the following equation:

where ${\stackrel{}{s}}_{k_{1q}}\left(t\right),{\stackrel{}{s}}_{k_{2q}}\left(t\right),...,{\stackrel{}{s}}_{k_{N_{s}q}}\left(t\right)$ represent the quantitative current density of ${\stackrel{}{s}}_{k_1}\left(t\right),{\stackrel{}{s}}_{k_2}\left(t\right),...,{\stackrel{}{s}}_{k_{N_s}}\left(t\right)$.

The following steps are performed on sLORETA-qm:

1. After the data were separately bandpass filtered in β-band or γ-band, the sLORETA-qm technique was performed for a certain time range. In this step, the lead field vector and matrix were calculated with a spherical model fit to the brain area of each subject. The center of the spherical model was determined as the mid point between anterior commissure and posterior commissure. The brain area in the spherical model was divided into approximately 2000 voxels with 10 × 10 × 10 mm³. The actual count of voxels depends on the head size of each subject. MEG data from all sensors were used for this calculation, thus avoiding arbitrary sensor selection.
2. The sLORETA-qm technique was performed on the results of the previous step. That is, voxels with spatially or locally maximum at each time point were selected and transformed using Eq. (8) to obtain quantitative values. In this step, voxels in areas without effect from brain activities, such as in the vicinity of eyes, were excluded. This operation can be performed automatically by referring to the spherical model of each subject.
3. Finally, calculated current density [nAm/mm²] for each voxel was averaged to the selected time interval and images of time-averaged sources were super-imposed on individual images from MRI.

Extension of sLORETA-qm Technique to Induced Oscillatory Changes

Next, the sLORETA-qm technique was extended to the induced oscillatory changes. Induced activities were considered as the current intensity variation in two states, i.e., baseline (resting) and active (task-related) states. From the results of sensor time-frequency analysis, the active state was determined to be the time range of β-band ERD or γ-band ERS appearing after the trigger, while the baseline state was determined to be that before the trigger with the same time range in each subject. With the sLORETA-qm technique, “Baseline Source Image” and “Active Source Image” were obtained in each trial (Fig. 2). Next, the “Differential Source Image” was obtained by subtracting the current intensity of the baseline source image from that of the active source image. This “Differential Source Image” was considered to be a source image of the oscillatory changes in each trial. In this process, the source image of oscillatory changes in each trial as well as the quantities in each voxel could be obtained. Next, “Differential Source Images” from all trials were averaged, and this image was designated as the “Original Image” (Fig. 3a).

Statistical Evaluation of the Spatial Filter Output

As the original source image obtained in the previous step could have contained unexpected artifacts and less repeatable components, statistical evaluation of the spatial filter output was needed. In the present study, a random permutation test with the null hypothesis was applied for statistical evaluation [10, 11].

All baseline and active images were permuted randomly, and “Permuted Baseline Source Image” and “Permuted Active Source Images” were obtained. By subtracting them and averaging, a “Permutation Image” with quantities in each voxel was obtained (Fig. 3b). This permutation process was repeated 1000 times and the distribution of quantities in these 1000 “Permutation Images” was obtained in each voxel. By comparing the quantities in the “Original image” with that in permutation distribution, the one-tailed probability P-values were obtained in each voxel. One example of permutation test for β-band ERD in the voxel near the contralateral primary motor cortex (M1) in subject 1 was demonstrated (Fig. 3c). In this voxel, the quantity was -3.62×10^-3 [nAm/mm²] and the uncorrected P-value was 0.016.

Next, multiple comparison correction was conducted to identify the activated voxel at the whole-brain level. In the present study, a procedure was adopted based on the false discovery rate (FDR) [12]. In the voxel noted above, the corrected P-value was 0.038. The significance level was defined as corrected P < 0.05 and ERD in this voxel was considered to be significant.

Using these two statistical steps, permutation test and correction of multiple comparisons with FDR, a statistically significant source image was obtained.