Graph fractional-order total variation EEG source reconstruction

EEG source imaging is able to reconstruct sources in the brain from scalp measurements with high temporal resolution. Due to the limited number of sensors, it is very challenging to locate the source accurately with high spatial resolution. Recently, several total variation (TV) based methods have been proposed to explore sparsity of the source spatial gradients, which is based on the assumption that the source is constant at each subregion. However, since the sources have more complex structures in practice, these methods have difficulty in recovering the current density variation and locating source peaks. To overcome this limitation, we propose a graph Fractional-Order Total Variation (gFOTV) based method, which provides the freedom to choose the smoothness order by imposing sparsity of the spatial fractional derivatives so that it locates source peaks accurately. The performance of gFOTV and various state-of-the-art methods is compared using a large amount of simulations and evaluated with several quantitative criteria. The results demonstrate the superior performance of gFOTV not only in spatial resolution but also in localization accuracy and total reconstruction accuracy.

Abstract-EEG source imaging is able to reconstruct sources in the brain from scalp measurements with high temporal resolution. Due to the limited number of sensors, it is very challenging to locate the source accurately with high spatial resolution. Recently, several total variation (TV) based methods have been proposed to explore sparsity of the source spatial gradients, which is based on the assumption that the source is constant at each subregion. However, since the sources have more complex structures in practice, these methods have difficulty in recovering the current density variation and locating source peaks. To overcome this limitation, we propose a graph Fractional-Order Total Variation (gFOTV) based method, which provides the freedom to choose the smoothness order by imposing sparsity of the spatial fractional derivatives so that it locates source peaks accurately. The performance of gFOTV and various state-of-the-art methods is compared using a large amount of simulations and evaluated with several quantitative criteria. The results demonstrate the superior performance of gFOTV not only in spatial resolution but also in localization accuracy and total reconstruction accuracy.
Index Terms-EEG source localization, inverse problem, fractional-order total variation, alternating direction method of multipliers

I. INTRODUCTION
Electroencephalography (EEG) is a valuable tool for studying the brain activity with excellent temporal resolution (∼ms) as well as supporting free body movement. EEG source imaging aims at localizing brain sources from scalp measurements by solving the corresponding inverse problem. However, due to the limited number of measurements (∼ 10 2 ) comparing with the large number of possible dipole source localizations (∼ 10 4 ), the EEG inverse problem is highly under-determined and has infinite solutions. By imposing various types of regularizations, the uniqueness of the solution can be guaranteed. So far, a large number of EEG reconstruction methods have been proposed to improve reconstruction resolution and localization accuracy by employing various regularization techniques. One of the most widely used methods is based on the Tikhonov regularization using the 2 -norm, such as minimum norm estimate (MNE) [1] and standardized low resolution brain electromagnetic tomography (sLORETA) [2]. The major drawback of *This work was supported in part by the California Capital Equity LLC and the Keck foundation. 1 Ying Li is with the Department of Bioengineering, University of California, Los Angeles, CA 90095, USA yingli.ucla@gmail.com this type of methods is its relatively low spatial resolution.
To further improve the spatial resolution, the 1 -norm based methods are proposed with sparsity constraint in the original source domain, such as minimum current estimate (MCE) [3]. Although the spatial resolution is greatly improved, these methods only involve a few elements in each source region, thus failing to identify the spatial extent of distributed sources. Recently, the total variation (TV) based methods [4], [5] are developed by imposing sparsity of the first spatial derivative rather than the original source. It shows great potential in preserving sharpness of edges and estimating the spatial extent of distributed sources. However, relying on the assumption that the underlying source is constant in each subregion, the TV-based methods suffer from the socalled staircase artifacts for the sources with higher order smoothness, which results in loss of the gradual intensity change in the reconstructed source and poor localization of source peaks.
In this paper, we propose a novel EEG source reconstruction method based on the Fractional-Order Total Variation (FOTV) to improve the reconstruction accuracy of the brain image. FOTV has been recently proposed to solve image processing problems [6], [7]. Different from TV which imposes sparsity of first-order spatial derivatives, FOTV can choose a more elegant smoothness order for the underlying source by imposing sparsity of α-order derivatives (α > 0). Therefore, the proposed method is capable of reconstructing the brain image with higher-order smoothness and preserving natural intensity changes of the brain image. As a consequence, the localization accuracy of peaks is greatly improved compared to TV-based methods. In order to extend the traditional FOTV defined on a 2D rectangular grid to an irregular cortex surface, we treat the cortex surface as a graph and define a novel graph FOTV (gFOTV) using shortest paths on the graph. In fact, TV-based methods can be considered as a special case in our proposed framework when α = 1. In addition, we derive an efficient algorithm using the alternating direction method of multipliers (ADMM). Finally, a large variety of simulations are conducted under different source sizes with a realistic head model. The proposed method is compared with several state-of-the-art ones both qualitatively and quantitatively, which demonstrates superior performance in terms of total reconstruction accuracy, spatial resolution and localization accuracy.

II. EEG INVERSE PROBLEM
Given acquired electrical recordings, the "EEG inverse problem" aims to identify the distributions of currents in a human brain. More precisely, after discretizing the cortex surface into M triangles (called "voxels"), the electrical potential b ∈ R N on the scalp measured by electrodes is related to the current density u where n denotes the noise, and the lead field matrix A can be calculated by constructing a head model using the geometric and electrical properties of the head. In this work we use the realistic head model and source model templates provided by Fieldtrip toolbox [8]. In general, since the number of voxels M (∼ 10 4 ) is much larger than that of electrodes N (∼ 10 2 ), the corresponding inverse problem is ill-posed. In order to obtain a unique solution for the problem, it is necessary to consider an appropriate regularization based on geometry of the cortex surface as well as the current density distribution of sources.
III. GRAPH FRACTIONAL-ORDER TOTAL VARIATION Fractional-order TV has been proposed recently in image processing to enhance smoothness with relatively low computational cost [6], [7], [9]. In this paper, we focus on the anisotropic version in which the objective function of the derived minimization problem is separable. Recall that the anisotropic fractional-order TV in a 2D rectangular mesh is defined as follows: where α ∈ (1, 2). Here the fractional derivative is based on the Grüwald-Letnikov derivative definition [10] (D α where the coefficients are It is easy to see that w α (0) = 1 and w α (1) = −α. Moreover, if α = 1, then TV α is the traditional TV. If α = 2, D α x /D α y approximates the second partial derivative of u along the x/y-direction. Although the above definition is also valid for α ∈ (0, 1)∪(2, ∞), it has experimentally shown that α ∈ (1, 2) achieves the best performance in applications [7]. Since the cortex surface is an irregular 3D surface consisting of gyrus and sulcus, we define a graph α-order TV with α ∈ [1, 2] tailored to such surface. After discretizing the cortex surface as mentioned in Section II, we create a graph whose nodes correspond to the centroids of all triangles. For a specific node v i , let d(v i , v j ) be the number of nodes on the shortest path connecting the nodes v i and v j , which is in or close to a geodesic of the underlying surface passing through v i and v j . Given a path p = (v i=m 0 , v m 1 , . . . , v m K ) where the shortest distance between the nodes v m 0 and v m j is j nodes, the fractional-order derivative along the path p is defined as Then the discretized fractional-order TV of u is defined as follows: where P(i; K) is the set of all paths starting from the ith node with length of K nodes. The shortest path between each node pair can be calculated by the breadth-first search (BFS) algorithm. For a specific node v i , the nodes at level k, i.e., the nodes have shortest distance k from v i , are assigned the weight w α (k). In particular, if α = 1, then D α is exactly the finite difference operator used in the TV-based methods. If α = 2, then w α (k) = 0 for k > 2 which implies that all nodes at level more than two are assigned zero weight. If α ∈ (1, 2), then the weights w α (k) will gradually decay as k goes to infinity. As the value of α increases from 1 to 2, the decay rate of weights w α (k) becomes larger. Note that K specifies the maximal level of nodes to be used. Due to the sparse structure of the underlying u, it is sufficient to choose K ≤ 4 levels of neighboring nodes to achieve the desired accuracy in our experiments.

IV. MODEL AND ALGORITHM
After defining FOTV on triangular mesh of cortex surface, we propose the following fractional-order TV regularized EEG source reconstruction model to improve high order smoothness of the brain image min u 1 2 Au where λ > 0 is regularization parameter, which controls the tradeoff between the data fidelity term and the sparsity term. By change of variables, the above problem can be converted to a linear equality constrained minimization problem min u,v Then the ADMM yields the following algorithm Here the parameters ρ > 0, γ ∈ (0, ( √ 5 + 1)/2) and the shrinkage operator is defined componentwise as shrink(u, µ) i = sign(u i ) max{|u i | − µ, 0}.

A. Simulation Protocol
In our simulation, Gaussian patches are used to simulate sources in the brain. To represent sources at different locations, we randomly select three sources located at different lobes of the cortex surface. In addition, to evaluate the performance of the proposed methods for different source sizes, we simulate sources containing 50, 100, 150, 200, 250 voxels, corresponding to 1.0cm, 1.4cm, 1.7cm, 1.9cm, 2.1cm in radius, respectively. After sources are generated, random independent and identically distributed (i.i.d.) Gaussian noise is added to each voxel as background neural noise. Additionally, electrode and electronic noise are added to each channel. For the signal-to-noise ratio (SNR), it is set to be 10dB by default. The number of electrode used is 346, and the number of voxels in the source model is 10240. Finally, the signal is normalized to between 10µV to 100µV , which is a typical range of amplitude for an adult EEG signal [12]. To reduce round-off errors, both the lead field matrix A and the electrical potential b are scaled by 10 5 .

B. Quantitative Criteria
The following three quantitative criteria are employed to evaluate performance of source reconstruction methods.
• Total reconstruction error (TRE), aiming to estimate the total reconstruction accuracy. It measures the relative error between the reconstructed source and the true one, and is defined as T RE = û − u 2 / u 2 , where u andû are the true and reconstructed source respectively. • Localization error (LE), aiming at evaluating the localization accuracy. It measures the averaged localization error between the peaks of the true and reconstructed sources [13], [14], and can be expressed as LE = ∑ k LE k /K, LE k = {d ki |i = argmin i ∈I k u i 2 },where d ki is the distance between the ith voxel to the peak of the kth true source, and I k is a set of the voxel indices that are closest to the peak of the kth source. • Degree of focalization (DF) is used to evaluate spatial resolution by measuring how focal the reconstructed source is. It is defined as DF = û S 2 2 / u S 2 2 [14], where u S is u restricted on the original patch area S.

C. Parameter Selection
In the operator D α , the parameter α specifies the order of the spatial derivative domain where we want to impose sparsity constraint on. As α becomes larger, we impose sparsity of higher order derivatives, and therefore the reconstructed sources become smoother and decay faster. When α = 1, the source decay speed is the slowest, which is the case of TV that the current density is piecewise constant. In our experience, it works well when α = 2. By choosing α to be a fraction between 1 and 2, all level nodes will be assigned a non-zero weight, which enhances reconstruction smoothness. Specifically, α = 1.6 is appropriate for all our experiments. In addition, the spatial resolution will be higher than that of α = 2. In the proposed Algorithm (4), the regularization parameter λ controls the balance between the data fidelity term and the sparsity term. When the source size becomes larger, λ should be tuned to be a little larger. As the noise level becomes higher, λ needs to be a little smaller. According to our experience, it works pretty well by simply choosing λ around 1. The parameter ρ affects the convergence speed of Algorithm (4), and is set as 10λ by default. Finally, the parameter γ is fixed as 1.

D. Simulation Results
First, we evaluate the performance of the proposed method using α = 2 and α = 1.6, and compare with several state-ofthe-art methods: sLORETA, 1 -norm based and TV-based methods. Fig.1 shows the reconstruction results of three sources at different locations of the cortex surface. One can see that the spatial resolution of sLORETA is very low, since it is based on the Tikhonov regularization using 2norm. The 1 -norm based method MCE greatly improves the spatial resolution because it imposes sparsity of the source in itself, but it yields over-focused reconstructed sources and fails to identify the spatial extent of the sources. The TVbased method shows good performance in identifying the spatial extent and preserving edges. However, it does not recover the varying intensity of the source, because of the assumption that the source intensity is piecewise constant. In contrast, the proposed method generates the reconstructed images closest to the ground truth. It not only provides high spatial resolution, but also successfully reconstructs the intensity variation and the spatial extent of the sources. It is worth noting that α = 1.6 is able to further enhance the spatial resolution compared to α = 2.
We further evaluate the proposed method quantitatively by using various quantitative criteria defined above. To show the robustness and consistency, we conduct a large number of simulation tests with different source sizes, and average results of 50 realizations with different noise configurations, illustrated in Fig.2. In the TRE curve, one can see that the proposed method provides the smallest total reconstruction error. Notice that the TV-based method also shows a relatively small total reconstruction error compared to other methods. In the LE curve, our proposed method shows the smallest localization error, since it makes the brain image smoother by imposing sparsity in higher order derivative instead of the first order derivative. As is well known that sLORETA can provide high localization accuracy, sLORETA also shows small localization error. TVbased method, however, produces relatively large localization error, which is because the reconstructed sources are constant in each subregion therefore it has difficulty pinpointing the peak of the source accurately. Finally, the DF curve shows that all of the methods with sparsity constraints show high focalization degree. It is interesting to notice that for gFOTV-based methods, α = 1.6 further enhances the spatial resolution. Compared to other methods, the spatial resolution of sLORETA is the lowest, as indicated by its name "standardized low resolution brain electromagnetic tomography". To sum up, our proposed method demonstrates superior performance from all the criteria including total reconstruction accuracy, localization accuracy and spatial resolution.

VI. CONCLUSIONS
In this paper, we propose an efficient and accurate EEG source reconstruction method by defining a novel graph fractional-order total variation (gFOTV) adapted for a triangular mesh of the cortical surface. This method imposes sparsity in α-order spatial derivatives with α ∈ [1, 2], which includes the TV-based methods as a special case when α = 1. By tuning the parameter α, the proposed method provides the freedom to choose a more elegant order of smoothness for the underlying brain image. Therefore, it not only provides high spatial resolution, but also recovers the current density variation and localizes the source peaks with high accuracy. In addition, the proposed algorithm is parameter friendly in the sense that the parameter selection is not sensitive to source size and noise level. A variety of simulation experiments demonstrate that the proposed method consistently outperforms the state-of-the-art methods both qualitatively and quantitatively. In future work, our proposed framework could be extended by using spatially variant smoothness orders for different subregions in an automatic way to further improve the reconstruction accuracy.