Coupling of Point Collocation Meshfree Method and FEM for EEG Forward Solver

Abstract

For solving electroencephalographic forward problem, coupled method of finite element method (FEM) and fast moving least square reproducing kernel method (FMLSRKM) which is a kind of meshfree method is proposed. Current source modeling for FEM is complicated, so source region is analyzed using meshfree method. First order of shape function is used for FEM and second order for FMLSRKM because FMLSRKM adopts point collocation scheme. Suggested method is tested using simple equation using 1-, 2-, and 3-dimensional models, and error tendency according to node distance is studied. In addition, electroencephalographic forward problem is solved using spherical head model. Proposed hybrid method can produce well-approximated solution.

1. Introduction

Electroencephalography (EEG) is a technique to measure and analyze the scalp potential [1] and widely used for clinical applications and scientific researches because EEG has noninvasiveness and high temporal resolution [2, 3]. As electric potential on scalp is a reflection of neuronal activation in a brain [4], solving inverse problem of EEG, or reconstructing distribution of current source, is an important mathematical tool to investigate brain functions and diseases [5]. Because precise solution of forward problem is essential for accurate inverse solution, we proposed a new forward solver with high accuracy in this paper.

Solving forward problem of EEG is to obtain potential distribution of head due to neuronal electric current [1, 6, 7]. Although there are many methods to solve forward problem, finite element method (FEM) is frequently used because it can consider realistic head shape, inhomogeneity, and anisotropy [7–9]. However, according to Awada’s research, variation of dipole location in a first-order element cannot affect forward solution, which causes poor spatial resolution of FEM [9]. Furthermore, location of dipole is also a problem in node-based distributed source model. To overcome these limitations, we adopted a meshfree method for region near dipole.

Because meshfree methods do not need finite elements, mesh generation process is not required. So, difficult situation to generate finite elements such as small objects inside a large analysis domain [10], crack propagation [11], and moving objects [12] can be analyzed by meshfree methods. However, meshfree methods consume longer time to solve a system matrix than FEM does, because the number of nodes related with a certain node of meshfree method is more than that of FEM.

In this paper, a new hybrid method of meshfree method and FEM for forward solver is proposed. There were some tries to surmount disadvantages of FEM and meshfree method by combining two methods in previous studies [13–16]. While most meshfree methods with weak formulation and integration cells are studied for coupling with FEM [13–16], point collocation scheme is applied for this study because of its efficiency [17, 18]. Among many sorts of meshfree methods [19–23], fast moving least square reproducing kernel method (FMLSRKM) is adopted because FMLSRKM can evaluate all approximated derivatives of shape functions up to the order of the basis polynomials simultaneously, without additional computations [16, 17]. We explain FMLSRKM and its coupling with FEM in the next section. Then, using simple models, the suggested method is tested and error is evaluated. Moreover, EEG problem with 2D model is treated to show that the suggested method can be a forward solver with small error.

2. Combination of FEM and FMLSRKM

2.1. FMLSRKM

In FMLSRKM [10, 16, 17], an approximated solution 𝑢^ℎ at 𝐱 near 𝐱 can be determined using polynomial bases 𝐏_𝑚 and their coefficient vector 𝐚(𝐱), like𝑢^ℎ𝐱,𝐱=𝐏_𝑚𝐱−𝐱𝜌⋅𝐚𝐱,(2.1) where 𝑚 is the highest order of polynomial and 𝜌 is the dilation parameter which represents region of influence of 𝑢^ℎ(𝐱,𝐱). In order to obtain coefficient vector, localized error residual functional𝐽𝐚𝐱=^NP_𝐼=1||||𝑢𝐱_𝐼−𝐏_𝑚𝐱−𝐱𝜌⋅𝐚𝐱||||²𝜌1^𝑛Φ𝐱^𝐼−𝐱𝜌,(2.2) should be minimized. Hence,𝐚𝐱=𝑀⁻¹𝐱^NP_𝐼=1𝐏_𝑚𝐱−𝐱𝜌𝑢𝐱^𝐼Φ_𝜌𝐱_𝐼𝑀−𝐱,(2.3)𝐱=^NP_𝐼=1𝐏_𝑚𝐱−𝐱𝜌𝐏^𝑇𝑚𝐱−𝐱𝜌Φ^𝜌𝐱_𝐼−𝐱,(2.4) where NP is the number of nodes in the local area. So, (2.1) can be rewritten as𝑢^ℎ𝐱,𝐱=𝐏^𝑇𝑚𝐱−𝐱𝜌𝑀⁻¹𝐱^NP_𝐼=1𝐏_𝑚𝐱−𝐱𝜌Φ^𝜌𝐱_𝐼−𝐱𝑢𝐱_𝐼,(2.5) whereΦ_𝜌𝐱_𝐼−𝐱=1𝜌^𝑛Φ𝐱^𝐼−𝐱𝜌(2.6) is called window function, whereΦ(𝐱)=(1−‖𝐱‖)^𝑡0,,wothheenrw‖is𝐱e‖.<1,𝑡>0,(2.7)

After moving least square scheme is applied, approximated solution (2.1) can be rewritten as𝑢(𝐱)=^NP_𝐼=1Ψ^𝐼[𝛼]𝑢𝐱_𝐼,(2.8) whereΨ^𝐼[𝛼]=𝛼!𝜌^|𝛼|𝐞^𝑇𝛼𝑀⁻¹𝐱𝐏_𝑚𝐱−𝐱𝜌Φ^𝜌𝐱_𝐼−𝐱,(2.9) is the shape function of FMLSRKM.

2.2. Coupling of Shape Functions

For the sake of an explanation, in this section, 1-dimensional shape functions are considered. Because point collocation method is adopted for FMLSRKM, quadratic shape function is used. In FEM, influence of a shape function is confined in an element (see Figure 1(a)). In addition, the value of shape function is 1 at the node and 0 at the other node, and the sum of shape functions at any point is 1, which makes shape function the partition of unity. However, in the case of FMLSRKM, an area covered by a shape function of a node is usually larger than that of FEM (see Figure 1(b)). Furthermore, the value of shape function at the node is not 1 and shape function is not the partition of unity, so the approximated value inside two nodes is interpolated by more than two shape functions in many cases.

(a)

(b)

Figure 1 

Shape functions of (a) FEM and (b) FMLSRKM. A shape function of FEM is a partition of unity, but that of FMLSRKM is not.

So, an interfacial area between FEM domain and FMLSRKM domain has one finite element shape function and many FMLSRKM shape functions. Hence, adequate weighting functions should be considered for shape functions of the hybrid method. While linear weighting functions are applied on both FEM and meshfree in the previous study [15], step functions are used for weight functions as shown in Figures 2(a) and 2(b). Consequently, Figure 2(c) is shape functions for coupling method, because FMLSRKM with point collocation scheme uses quadratic polynomials as bases.

(a)

(b)

(c)

Figure 2 

Combination of shape functions in an interfacial area of FEM and FMLSRKM. (a) Shape functions of FEM, (b) shape functions of FMLSRKM, and (c) shape functions of hybrid method. Shape function of FMLSRKM is fully adopted as shape function of the interfacial area (inside two vertical dotted lines).

The suggested hybrid method has some advantages. First, the system matrix of the hybrid method is sparse almost same as that of FEM, if most of analysis domain can remaine as finite element and meshfree region is well defined and specific. Second, adaptive approach is easy for meshfree region. In many cases, meshfree region may be a high-error area, because a region which is difficult to analyze is selected as meshfree region. Third, point collocation scheme is used for FMLSRKM, which makes applying boundary condition easy. Forth, the region of FMLSRKM has less error, since FMLSRKM of this paper has quadratic shape function.

3. Test and Result

Tests are performed using 1-, 2-, and 3-dimensional models with∇²𝜙=−1.(3.1) 1D test model is shown in Figure 3(a). In this figure, a domain for FMLSRKM Ω_FM is inside FEM domain Ω_FEM, and their interfacial area Ω_IN is surrounded by Γ_FEM and Γ_FM. At the ends, Dirichlet boundary conditions are applied as 𝜙=0. In this test problem all boundaries are in FEM region. It is also possible that FMLSRKM region has boundaries, and applying boundary condition into FMLSRKM is even easier because FMLSRKM of this paper uses collocation scheme.

(a)

(b)

(c)

Figure 3 

Test problems. In these problems, numbers, variables, and solutions are unitless. (a) 1D, (b) 2D. The solid line along 𝑦=5 is the cross-sectional line for error evaluation, and (c) 3D.

In Figures 3(b) and 3(c), test problems in 2D and 3D are shown, respectively. In these problems, domains for FMLSRKM are in the middle of the whole domains, and their boundary conditions are applied as shown in Figures 3(b) and 3(c), respectively. Hence, it is expected that the solutions of Figures 3(a), 3(b), and 3(c) are the same. Furthermore, these problems have analytic solution as𝜙=−12𝑥²+10𝑥,(3.2) which can be compared to approximated solutions. The results of the hybrid method for Figures 3(a), 3(b), and 3(c) are shown in Figures 4(a), 4(b), and 4(c), respectively. All solutions are well approximated.

(a)

(b)

(c)

Figure 4 

The results of the hybrid method. (a), (b), and (c) are the results of Figures 3(a), 3(b), and 3(c), respectively.

4. Error Study

In this paper, error was evaluated using various intervals of nodes, or ℎ. For evaluation of errors, the solid line along 𝑦=5 of Figure 3(b) is used. An error is defined as a difference between the exact solution and an approximated solution by the hybrid method. In Figure 4(a), it seems that the hybrid method generates the same result as the exact solution (3.2). However, actual error is not zero as shown in Figure 5. FEM region has higher error than FMLSRKM region because FEM uses 1st-order polynomials as shape functions while FMLSRKM uses 2nd-order ones. In FMLSRKM region, the error is nearly zero because 2nd-order FMLSRKM can describe (2.2) which is 2nd-order polynomials.

Figure 5 

The error of approximation by the hybrid method. The error of FEM is higher than that of FMLSRKM because of the orders of bases for two methods.

To investigate the error of FMLSRKM region, inhomogeneous forcing term should be applied. So,∇²𝜙=−0.6𝑥,(4.1) is considered and the boundary conditions are the same as the previous example. In this case,𝜙=−0.1𝑥³+40𝑥,(4.2) is the exact solution of (2.3). For this study, defining ℎ as the interval of nodes, ℎ=1, ℎ=0.5, and ℎ=0.25 are used. The case of ℎ=1 is as in Figure 3(b), and Figure 6(a) represents its absolute error which is evaluated along 𝑦=5. In FMLSRKM region, the error at the middle of each interval is zero, while the error is not zero but the maximum in FEM region. Figure 7 explains the reason, and this phenomenon is caused by quadratic approximation of FMLSRKM. As shown in Figure 7(a), linear approximation of FEM is not capable of following the cubic solution and produces maximum error at the middle. On the other hand, because the quadratic approximation of FMLSRKM and the exact solution are the same at the two nodes (𝑥=7 and 8 for Figure 7(b)), the quadratic approximation crosses the exact solution, which makes the error zero at the middle. For the cases of ℎ=0.5 and 0.25, the absolute errors are shown in Figures 6(b) and 6(c). The errors become smaller as the nodes become closer to next nodes. From Figure 6, tendency of relative errors is shown in Figure 8. Axes are drawn in logarithmic scale, and Figure 8 shows that the hybrid method has good characteristic of convergence.

(a)

(b)

(c)

Figure 6 

Errors of the hybrid method along 𝑦=5 of Figure 3(b). (a) ℎ=1, (b) ℎ=0.5, and (c) ℎ=0.25. The values of 𝑦-axis are absolute errors which decrease according to the size of ℎ.

(a)

(b)

Figure 7 

Examples of FEM and FMLSRKM. (a) Linear approximation of FEM shows maximum difference at the middle of the interval, and (b) Quadratic approximation of FMLSRKM crosses the exact solution.

Figure 8 

Relative error of FEM and FMLSRKM according to the interval of nodes.

5. EEG Forward Problem

Governing equation and boundary condition of EEG forward problem are∇𝐧⋅(⋅𝜎𝜎∇∇𝑉𝑉)==0∇⋅o𝐉n_𝑝𝜕iΩn,Ω,(5.1) where 𝜎, 𝑉, 𝐉_𝑝, and 𝐧 are electric conductivity, potential, and primary current dipole, and unit vector normal to boundary, respectively. Test model and result are shown in Figure 9. Simple spherical head model as shown in Figure 9(a) is used, and a current dipole is in FMLSRKM region. Because divergence of primary current of (5.1) can be represented as combination of one source and one sink [9], a node at dipole location should be divided into 2 nodes. In Figures 9(b) and 9(c), brighter region has larger error. Those figures show that the hybrid method can approximate the forward solution better than FEM.

(a)

(b)

(c)

Figure 9 

Test for EEG forward problem. (a) A dipole shown as an arrow is in the middle of FMLSRKM region, (b) result using FEM, and (c) result using the suggested hybrid method.

6. Conclusion

A hybrid method of FMLSRKM and FEM is suggested in this paper, and it shows good result of 1-, 2-, and 3-dimensional boundary value problems. In addition, error of the hybrid method is studied, and proposed method shows good convergence in 2-dimensional model. In Figure 8, relative errors in amplitude are compared between FEM and FMLSRKM. From this figure, we want to say not that FMLSRKM is better method than FEM but that the hybrid method does not lose own error characteristic of each method.

We also tested the hybrid method for EEG forward problem and verified that the method can produce small-error solution. Actually, it is possible to reduce error of solution by hp-FEM [24]. However, adaptive strategy can be more adequate for meshfree methods than FEM, because new nodes are just placed on regions which have high error without generating new mesh structures [10, 19]. Hence, error estimation and adaptive scheme will be applied onto the proposed method.

This research shows that Poisson and Laplace equations were successfully handled by the suggested hybrid method. Therefore, realistic problems can be treated such as moving objects and singular points, and it is expected that the hybrid method will take its own role for problems with special conditions as well as EEG forward problem.

Acknowledgment

This study was supported by a grant from the Korea Healthcare Technology R and D Project, Ministry for Health, Welfare and Family Affairs, Republic of Korea (A090794).