A New Weighted Correlation Coefficient Method to Evaluate Reconstructed Brain Electrical Sources

Abstract

Various inverse algorithms have been proposed to estimate brain electrical activities with magnetoencephalography (MEG) and electroencephalography (EEG). To validate and compare the performances of inverse algorithms, many researchers have used artificially constructed EEG and MEG datasets. When the artificial sources are reconstructed on the cortical surface, accuracy of the source estimates has been difficult to evaluate. In this paper, we suggest a new measure to evaluate the reconstructed EEG/MEG cortical sources more accurately. To validate the usefulness of the proposed method, comparison between conventional and proposed evaluation metrics was conducted using artificial cortical sources simulated under different noise conditions. The simulation results demonstrated that only the proposed method could reflect the source space geometry regardless of the number of source peaks.

1. Introduction

Noninvasive measurements of brain electrical activities with electroencephalography (EEG) and magnetoencephalography (MEG) enabled us to estimate the underlying cortical activities, thereby contributing to the rapid development of clinical and cognitive neuroscience. To estimate the cortical electrical activities from EEG and MEG, of which the process is often called EEG/MEG source imaging, highly underdetermined inverse problems have to be solved using linear or nonlinear inverse algorithms since the source estimation from EEG and MEG signals is an ill-posed problem, which generally produces blurry or inaccurately positioned source estimates [1]. Many mathematical approaches and techniques have been proposed to estimate accurate source locations and strengths. Among them, minimum-norm estimate (MNE) has been the most widely studied inverse algorithm as MNE is simple and has linearity [2]. MNE chooses a source distribution where the 𝑙2 norm of the current distribution is minimized. On the contrary, minimum current estimate (MCE) selects a source where the 𝑙1 norm of the current is minimized [3]. Other than those two representative algorithms, there have been several modifications of norm minimization, for example, low-resolution electrical tomography (LORETA) algorithm [4] and the focal underdetermined system solution (FOCUSS) algorithm [5].

When a new source imaging algorithm is proposed, the performance of the inverse algorithm needs to be verified and compared with those of the existing ones. Since the reconstructed source distributions are hard to be verified using in vivo experiments, many researchers have used artificial EEG/MEG human skull phantoms [6] or realistically simulated EEG/MEG datasets. Since the use of simulated EEG/MEG data allows us to readily adjust and control noise levels and source configurations, that is, the number and size of source patches, most inverse algorithms are generally verified using simulated EEG/MEG data [2–5]. In recent simulation studies, the source spaces are generally constrained only on the interface between white and gray matter of the cerebral cortex, generally called cortical surface, considering neurophysiology. The orientations of the cortical sources are also assumed to be perpendicular to the cortical surface [7]. In such simulations, both the original source patches and the reconstructed sources are commonly distributed on the cortical surface generally tessellated with surface triangular elements.

For the evaluation of the reconstructed sources, evaluation metrics or error metrics need to be introduced to probe the similarity between the simulated and reconstructed sources. The well-known evaluation metrics are root mean square error (RMSE), shift of the maximum (𝑆_max), shift of the center of mass (𝑆_cm), and the correlation coefficient (CC) [6]. Each metric has its own pros and cons. In contrast to the conventional geometric error metrics such as 𝑆_max and 𝑆_cm, RMSE and CC do not reflect the geometry of the cortical surface. However, compared to 𝑆_max and 𝑆_cm, RMSE and CC are reliable specifically when the source distributions are not concentrated to a single peak. For more accurate and robust estimation of the accuracy of reconstructed EEG/MEG sources, we modified CC by giving the geodesic distance weights to the reconstructed sources to reflect the geometric information of cortical surface. To validate the new evaluation metric, named weighted correlation coefficient (WCC), some representative examples were used.

2. Methods

2.1. EEG/MEG Inverse Problem

When a set of 𝑛 possible source locations and 𝑚 sensor positions is given, thanks to the linearity of Maxwell’s equations, an EEG/MEG forward model can be described as 𝐛=𝐊𝐣+𝐬, where 𝐊 is an 𝑚 by 𝑛 EEG/MEG lead field matrix, 𝐣 is an 𝑛 by 1 unknown source vector, 𝐛 is an 𝑚 by 1 recorded EEG/MEG data, and s is the additional sensor noises. The inverse problem for estimating 𝐣 from 𝐛 has no unique solution. To estimate the possible solutions, MNE adopts the following minimization problem: min_j‖𝐣‖₂subjectto𝐛=𝐊𝐣+𝐬. (2.1)

Then, the estimated solution 𝐣̃ can be written as 𝐣̃=𝐊^𝑇𝐊𝐊^𝑇+𝜆𝐈⁻¹𝐛, (2.2) where 𝜆 is a regularization parameter, which was determined using the generalized cross-validation method [8].

2.2. Conventional Evaluation Metrics

We assume that both the simulated true sources j and the estimated sources j̃ are distributed on the 3D cortical surface. The dimension of both vectors is 𝑛 by 1, where 𝑛 is the number of nodes on the source space. We firstly summarize four conventional evaluation metrics, having been frequently used for assessing the accuracy of the source estimates.

2.2.1. Root Mean Square Error

The root mean square error (RMSE) is the most well-known and convenient way to measure the error between the actual source and the estimated source. RMSE is formulated as RMSE=1𝑛^𝑛_𝑖=1𝑗_𝑖−𝑗̃_𝑖², (2.3) where 𝑗_𝑖 and j_ĩ are the 𝑖th elements of j andj̃respectively.

This metric is easy to implement and can be used regardless of the shapes of the source distributions. However, RMSE does not reflect the geometry of the cortical surface since RMSE is computed with just vectored values.

2.2.2. Shift of the Maximum of the Estimate

The shift of the maximum of the estimate (𝑆_max) is the simplest measure which reflects the geometry of the source space. 𝑆_max indicates the distance between the locations where the maximum intensities of sources are generated. The maximum intensities of the actual and reconstructed source are assumed to be located at 𝐫_𝐦𝐚𝐱 and 𝐫̃_𝐦𝐚𝐱; respectively, 𝐫_𝐦𝐚𝐱=max_𝐫𝑖𝑗_𝑖,𝐫̃_𝐦𝐚𝐱=max_𝐫𝑖𝑗_𝑖̃, (2.4) where 𝐫_𝑖 is the coordinate of 𝑖th node, then 𝑆_max is defined as 𝑆_max=𝐫_𝐦𝐚𝐱−𝐫̃_𝐦𝐚𝐱² (2.5) and ranged from 0 to 𝑑_max, the maximum distance within the brain.

This measure is reliable only when the actual source is concentrated around the location of the maximum source intensity because it does not consider the distributions of the cortical sources. When 𝑆_max is adopted as a measure, the merit of distributed source modeling disappears. For example, even when the extents of the true source and the reconstructed sources are largely different, identical maximum location makes the 𝑆_max value be 0.

2.2.3. Shift of the Center of Mass

The center of mass has been widely used for evaluating various algorithms adopted not only in EEG and MEG but also in other functional brain imaging techniques such as functional magnetic resonance imaging (fMRI) and positron emission tomography (PET). The center of mass of the actual source 𝐫_𝐜𝐦 and the center of mass of the reconstructed source 𝐫̃_𝐜𝐦 are computed as 𝐫_𝐜𝐦=∑^𝑛_𝑖=1𝑗_𝑖𝐫_𝑖∑^𝑛_𝑖=1𝑗_𝑖,𝐫̃_𝐜𝐦=∑^𝑛_𝑖=1𝑗_𝑖̃𝐫_𝑖∑^𝑛_𝑖=1𝑗_𝑖̃. (2.6)

As assuming the distributed source to be a dipole source placed on the center of mass of the source, the shift of center of mass (𝑆_cm) is defined as the distance between 𝐫_𝐜𝐦 and 𝐫̃_𝐜𝐦:𝑆_cm=𝐫_𝐜𝐦−𝐫̃_𝐜𝐦². (2.7)

𝑆_cm is similar to 𝑆_max in that the distributed source is considered as a point source placed at a single location. Therefore, 𝑆_cm is also reliable only when the simulated source is concentrated around 𝐫_𝐜𝐦. If the distribution of the source has a radial symmetry, 𝑆_cm becomes equivalent to 𝑆_max.

2.2.4. Correlation Coefficient

The correlation coefficient (CC), a concept adopted from statistics, is a measure of linear dependency between two variables, and the value ranges between −1 and 1. It has been widely employed as a standard measure in various fields of engineering and sciences. The conventional CC is defined as the covariance of 𝐣 and 𝐣̃ divided by the product of their standard deviations: 𝐣̃CC=cov𝐣,cov(𝐣,𝐣)cov𝐣,̃𝐣̃, (2.8) where the covariance is defined as cov𝐣,𝐣̃=1𝑛^𝑛_𝑖=1𝑗_𝑖−𝑗^∗𝑗_𝑖̃−𝑗^∗̃, (2.9) and 𝑗^∗ represents the mean value of the source j: 𝑗^∗=1𝑛^𝑛_𝑖=1𝑗_𝑖,𝑗^∗̃=1𝑛^𝑛_𝑖=1𝑗_𝑖̃. (2.10)

If the distribution of the reconstructed sources is similar to that of the actual sources, the value of CC is close to 1; if the distribution of the reconstructed sources is different from that of the actual sources, CC is close to −1. CC is reliable even when the source distribution is not concentrated to a single location or when the true source has many distinct peaks. However, similar to RMSE, CC cannot reflect the real geometry of the cortical surface.

2.3. New Algorithm: Weighted Correlation Coefficient

When we categorize four conventional measures mentioned in Section 2.2 in terms of geometric consideration, contrast to 𝑆_max and 𝑆_cm, RMSE and CC do not reflect the geometry of the cortical surface. However, RMSE and CC are reliable compared to 𝑆_max and 𝑆_cm when the source distribution is not concentrated to a single peak. To combine the advantages of both types of conventional measures, we modified CC by giving the source vector a weight reflecting geometrical information of cortical surface. The new evaluation measure, named weighted correlation coefficient (WCC), is defined as 𝐣̃WCC=cov𝐖𝐣,𝐖cov(𝐖𝐣,𝐖𝐣)cov𝐖𝐣,̃𝐖𝐣̃, (2.11) and 𝐖 is an 𝑛 by 𝑛 weighting matrix that can be computed as𝑑𝐖=^max𝐈_𝐧−𝐃𝑑_max, (2.12) where 𝐈_𝐧 is an 𝑛 by 𝑛 identity matrix. 𝐃 is an 𝑛 by 𝑛 distance matrix whose element is given as 𝐷_𝑖𝑗=𝐫_𝑖−𝐫_𝑗𝑘 (2.13) and 𝑑_max is the maximum value in 𝐃. If 𝑘=2, the Euclidean distance is employed, and if 𝑘=geo, then the geodesic distance is employed to obtain the distance matrix. The geodesic distance was computed by solving the Eikonal equation on the tessellated cortical surface [9]. The main diagonal of the weight matrix 𝐖 was filled with 1, and the off-diagonal elements were filled with values between 0 and 1. By multiplying weight matrix 𝐖 to the source vector j, the geometric information of cortical surface is considered.

Additionally, Euclidean or geodesic distance can be employed in the definition of the distance matrix 𝐃. Since the cortical surface of a human brain is folded, the geodesic distance is more suitable to reflect the geometric information of the cortical surface than the Euclidean distance. The Euclidian distance is computed by the Cartesian coordinates regardless of the geometrical feature of the cortical surface. However, as the geodesic distance implies the minimum distance along the surface, the geodesic distance between the two adjacent gyri should be greater than the Euclidian distance. Figure 1 is an example of the Euclidean and geodesic distance between each cortical surface vertex and a reference point located at right dorsolateral prefrontal cortex, corresponding to a column of the distance matrix.

(a)

(b)

Figure 1 

One column of (a) Euclidean and (b) geodesic distance matrix visualized on the cortical surface.

Once the distance matrix is computed, the weighting matrix is determined from (2.13). Figure 2 shows the Euclidean and geodesic weights corresponding to the Euclidean and geodesic distance matrices exampled in Figure 1. In contrast to Figure 1, as getting farther from the reference point, the weight is getting smaller.

(a)

(b)

Figure 2 

One column of (a) Euclidean and (b) geodesic weighted matrix visualized on the cortical surface.

The characteristics of the conventional and proposed measures are summarized in Table 1. Low values of RMSE, 𝑆_max, or 𝑆_cm and high values of CC and WCC indicate the accurate reconstruction. Only WCC is applicable to the case of multipeak and can consider the geometry of source space.

Table 1 

The merits and demerits of measures.

3. Results

To compare and verify the conventional and proposed measures, a simple two-dimensional example was simulated as shown in Figure 3. The source space was defined as a two-dimensional rectangle. The actual source distribution 𝐱 is given in Figure 3(a), and five reconstructed sources are given in Figures 3(b)–3(f), each of which was denoted as 𝐲_𝟏, 𝐲_𝟐, 𝐲_𝟑, 𝐲_𝟒, and 𝐲_𝟓. The source current intensities are indicated with different colors. If we evaluate the reconstructed sources based on visual inspection, anyone would agree that 𝐲_𝟏 is the most accurate reconstruction and 𝐲_𝟐 is the second best one. 𝐲_𝟓 seem to be the worst reconstruction as the peak location is farthest from the actual one and no reconstructed source is overlapped with the actual one. 𝐲_𝟑 and 𝐲_𝟒 seems to be better matched than 𝐲_𝟓, but it is difficult to judge which result is better. The result 𝐲_𝟑 has no commonly activated region with the actual source, but the distribution is close to the actual source distribution, whereas 𝐲_𝟒 has slightly overlapped region, but other regions are located far from the actual source location. If we assume visual inspection (VI) as a qualitative measure, the rank of the reconstructed sources can be expressed as VI(𝐱,𝐲_𝟏)>VI(𝐱,𝐲_𝟐)>VI(𝐱,𝐲_𝟑)≥VI(𝐱,𝐲_𝟒)>VI(𝐱,𝐲_𝟓).

(a) x

(b) 𝐲𝟏

(c) 𝐲𝟐

(d) 𝐲𝟑

(e) 𝐲𝟒

(f) 𝐲𝟓

Figure 3 

Example of a simulated two-dimensional source space: actual source distribution (a) and reconstructed sources (b)–(f).

We then employed the conventional and proposed quantitative measures for the evaluation of the reconstructions depicted in Figure 1 and summarized the result in Table 2. All measures commonly indicated that 𝐲_𝟏 is the best reconstruction and 𝐲_𝟓 is the worst reconstruction. However, the different metrics showed different evaluation results for 𝐲_𝟐, 𝐲_𝟑, and 𝐲_𝟒. In the case of RMSE, 𝐲_𝟑 was evaluated as the worst reconstruction, and 𝐲_𝟒 and 𝐲_𝟐 had an identical RMSE value, which was because RMSE was affected by the commonly activated regions regardless of the source geometry. In the case of 𝑆_max, which considers only the maximum location of the source, the results of 𝐲_𝟐 and 𝐲_𝟒 were equivalent. Similar to RMSE, CC classified 𝐲_𝟑 as the worst reconstruction and 𝐲_𝟒 and 𝐲_𝟐 had an identical CC value. Both 𝑆_cm and WCC evaluated the reconstruction results identically to the visual inspection results. However, if the actual source has multiple peaks, 𝑆_cm cannot be accurately evaluated.

Table 2 

Evaluation of reconstructions depicted in Figure 1.

We then applied the conventional and proposed evaluation metrics to the evaluation of distributed sources on the cortical surface. The cortical surface was extracted from a structural MRI of a standard brain atlas provided by the Montreal Neurological Institute (MNI). To extract and tessellate the cortical surface, we used CURRY6 for Windows (Compumedics, Inc., El Paso, TX). The actual source 𝐱 defined on the cortical surface is given in Figure 4(a), and the reconstructed sources are given in Figures 4(b)–4(f), each of which was denoted as 𝐲_𝟏, 𝐲_𝟐, 𝐲_𝟑, 𝐲_𝟒 and 𝐲_𝟓. If we verify the reconstructed sources by visual inspection, 𝐲_𝟏 seems to be the most accurate reconstruction and 𝐲_𝟐 seems to be the second best one. The two distributions 𝐲_𝟒 and 𝐲_𝟓 are very different from the actual source distribution. We can roughly estimate the rank of the reconstructed sources as VI(𝐱,𝐲_𝟏)>VI(𝐱,𝐲_𝟐)>VI(𝐱,𝐲_𝟑)>VI(𝐱,𝐲_𝟒)≥VI(𝐱,𝐲_𝟓).

(a) x

(b) 𝐲𝟏

(c) 𝐲𝟐

(d) 𝐲𝟑

(e) 𝐲𝟒

(f) 𝐲𝟓

Figure 4 

Simulation of three-dimensional cortical sources: actual source (a) and reconstructed sources (b)–(f).

The quantitative evaluation results with conventional and proposed metrics are shown in Table 3. In the case of RMSE, the RMSE values corresponding to 𝐲_𝟏~𝐲_𝟒 were increasing, which coincided with the visual inspection results, but the increment was very small compared to the absolute values of RMSE. The results of 𝑆_max,𝐸 and 𝑆_max,𝐺 showed that 𝐲_𝟓 is the better reconstruction than 𝐲_𝟑. Moreover, since the center of mass of 𝐲_𝟓 was located near the actual source, 𝑆_cm,𝐸 and 𝑆_cm,𝐺 of 𝐲_𝟓 were less than those of 𝐲_𝟐. Subscripted 𝐸 and 𝐺 indicate that the Euclidean and geodesic distance matrices were adopted, respectively. CC also could not distinguish the difference between 𝐲_𝟑 and 𝐲_𝟒. WCC_𝐸 and WCC_𝐺 are the WCC results when the Euclidean and geodesic distance matrices were adopted, respectively, to construct the weighting matrix. Both WCC_𝐸 and WCC_𝐺 were proved to be reliable since both results coincided well with the visual inspection results. However, compared to WCC_𝐺, WCC_𝐸 could not reflect the large difference between 𝐲_𝟐 and 𝐲_𝟑, which are located even in different hemispheres.

Table 3 

Evaluation of reconstructions depicted in Figure 4.

We performed extensive computer simulations to quantitatively compare the performance of the proposed measure with that of conventional measures. 2,000 locations on the cortical surface are selected randomly, and on each location a single constant source patch is generated. Consequently, our computer simulations used 2,000 cortical patches whose averaged radius is 6 mm (±1.2 mm). After solving the MEG forward problem with each cortical patch, different-level white Gaussian noise are added to the simulated MEG signal data. We set the signal-to-noise ratio (SNR) values from −10 dB to 30 dB. The reconstructed source is computed by the minimum norm method with each MEG signal then evaluated with different measures. Table 4 shows the averaged accuracies with conventional and proposed measures with respect to SNR. We expect that results of reconstruction to become more accurate as SNR is getting higher. In the case of RMSE, though the RMSE is increasing as SNR becomes lower, the difference between 10, 20, and 30 dB cases is not clear as much the cases of 0 and −10 dB. In the cases of geodesic measures (𝑆_max,𝐸, 𝑆_cm,𝐸, 𝑆_max,𝐺, 𝑆_cm,𝐺), the relation of accuracy and SNR is not consistent with the expected tendency. Moreover, the accuracy of low SNR is occasionally greater than that of high SNR in geodesic measures. In the case of CC with high SNR (10, 0, and −10 dB), there is only a marginal difference compared to that of low SNR. Only the accuracy measured by WCC is consistently decreasing as the SNR becomes lower.

Table 4 

Evaluation of reconstructions with 2,000 randomly selected source patches.

4. Conclusion

The geometric measures (𝑆_max, 𝑆_cm) could reflect the geometry of the source space, while the statistical measures (RMSE, CC) could be applied regardless of the distribution characteristics of the sources. In this paper, a new evaluation metric named WCC was proposed to combine the advantages of both types of evaluation metrics and showed enhanced performances compared to the conventional metrics. From the extensive simulation, we could conclude that the proposed measure is very promising to evaluate accuracy of reconstructed sources or EEG/MEG inverse algorithms.