Journal of Applied MathematicsVolume 2012 (2012), Article ID 840593, 9 pagesdoi:10.1155/2012/840593

Research Article

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

Chany Lee,¹ Jong-Ho Choi,² Ki-Young Jung,¹ and Hyun-Kyo Jung¹

¹College of Medicine, Korea University, 126-1 Anam-dong, Sungbuk-gu, Seoul 136-705, Republic of Korea²School of Electrical Engineering and Computer Science, Seoul National University, Seoul 151-742, Republic of Korea

Received 15 November 2011; Accepted 28 December 2011

Academic Editor: Chang-Hwan Im

Copyright © 2012 Chany Lee et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 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 should be minimized. Hence, where NP is the number of nodes in the local area. So, (2.1) can be rewritten as where is called window function, where

After moving least square scheme is applied, approximated solution (2.1) can be rewritten as where 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 1D test model is shown in Figure 3(a). In this figure, a domain for FMLSRKM is inside FEM domain , and their interfacial area is surrounded by and . At the ends, Dirichlet boundary conditions are applied as . 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 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 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 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, is considered and the boundary conditions are the same as the previous example. In this case, is the exact solution of (2.3). For this study, defining as the interval of nodes, , , and are used. The case of is as in Figure 3(b), and Figure 6(a) represents its absolute error which is evaluated along . 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 ( 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 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 of Figure 3(b). (a) , (b) , and (c) . 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 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).

References

H. Hallez, B. Vanrumste, R. Grech et al., “Review on solving the forward problem in EEG source analysis,” Journal of NeuroEngineering and Rehabilitation, vol. 4, no. 1, article 46, 2007.

A. M. Dale and M. I. Sereno, “Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach,” Journal of Cognitive Neuroscience, vol. 5, no. 2, pp. 162–176, 1993.

S. Baillet, J. C. Mosher, and R. M. Leahy, “Electromagnetic brain mapping,” IEEE Signal Processing Magazine, vol. 18, no. 6, pp. 14–30, 2001.

D. Purves, “Principles of cognitive neuroscience,” Sinauer Associates, vol. 83, no. 3, 757 pages, 2008.

C. H. Im, C. Lee, K. O. An, H. K. Jung, K. Y. Jung, and S. Y. Lee, “Precise estimation of brain electrical sources using anatomically constrained area source (ACAS) localization,” IEEE Transactions on Magnetics, vol. 43, no. 4, pp. 1713–1716, 2007.

Z. Zhang, “A fast method to compute surface potentials generated by dipoles within multilayer anisotropic spheres,” Physics in Medicine and Biology, vol. 40, no. 3, article 01, pp. 335–349, 1995.

J. C. Mosher, R. M. Leahy, and P. S. Lewis, “EEG and MEG: forward solutions for inverse methods,” IEEE Transactions on Biomedical Engineering, vol. 46, no. 3, pp. 245–259, 1999.

C. H. Im, C. Lee, H. K. Jung, Y. H. Lee, and S. Kuriki, “Magnetoencephalography cortical source imaging using spherical mapping,” IEEE Transactions on Magnetics, vol. 41, no. 5, pp. 1984–1987, 2005.

K. A. Awada, D. R. Jackson, J. T. Williams, D. R. Wilton, S. B. Baumann, and A. C. Papanicolaou, “Computational aspects of finite element modeling in EEG source localization,” IEEE Transactions on Biomedical Engineering, vol. 44, no. 8, pp. 736–752, 1997.

C. Lee, C.-H. Im, H.-K. Jung, H.-K. Kim, and D. W. Kim, “A posteriori error estimation and adaptive node refinement for fast moving least square reproducing kernel (FMLSRK) method,” Computer Modeling in Engineering & Sciences, vol. 20, no. 1, pp. 35–41, 2007.

T. Belytschko, Y. Y. Lu, and L. Gu, “Crack propagation by element-free Galerkin methods,” Engineering Fracture Mechanics, vol. 51, no. 2, pp. 295–315, 1995.

T. Rabczuk, R. Gracie, J.-H. Song, and T. Belytschko, “Immersed particle method for fluid-structure interaction,” International Journal for Numerical Methods in Engineering, vol. 81, no. 1, pp. 48–71, 2010.

T. J. Barth and M. Field, Meshfree Methods for Partial Differential Equations II, Springer, 2003.

B. N. Rao and S. Rahman, “A coupled meshless-finite element method for fracture analysis of cracks,” International Journal of Pressure Vessels and Piping, vol. 78, no. 9, pp. 647–657, 2001.

T. Rabczuk, S. P. Xiao, and M. Sauer, “Coupling of mesh-free methods with finite elements: basic concepts and test results,” Communications in Numerical Methods in Engineering, vol. 22, no. 10, pp. 1031–1065, 2006.

D. W. Kim and H. K. Kim, “Point collocation method based on the FMLSRK approximation for electromagnetic field analysis,” IEEE Transactions on Magnetics, vol. 40, no. 2, pp. 1029–1032, 2004.

D. W. Kim and Y. Kim, “Point collocation methods using the fast moving least-square reproducing kernel approximation,” International Journal for Numerical Methods in Engineering, vol. 56, no. 10, pp. 1445–1464, 2003.

C. Lee, D. W. Kim, S. H. Park, H. K. Kim, C. H. Im, and H. K. Jung, “Point collocation mesh-free method using FMLSRKM for solving axisymmetric Laplace equation,” IEEE Transactions on Magnetics, vol. 44, no. 6, Article ID 4526930, pp. 1234–1237, 2008.

T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free Galerkin methods,” International Journal for Numerical Methods in Engineering, vol. 37, no. 2, pp. 229–256, 1994.

W.-K. Liu, S. Li, and T. Belytschko, “Moving least-square reproducing kernel methods. I. Methodology and convergence,” Computer Methods in Applied Mechanics and Engineering, vol. 143, no. 1-2, pp. 113–154, 1997.

R. A. Gingold and J. J. Monaghan, “Smoothed particle hydrodynamics: theory and application to non-spherical stars,” Monthly Notices of the Royal Astronomical Society, vol. 181, no. 2, pp. 375–389, 1977.

W. L. Nicomedes, R. C. Mesquita, and F. J.S. Moreira, “A meshless local Petrov-Galerkin method for three-dimensional scalar problems,” IEEE Transactions on Magnetics, vol. 47, no. 5, pp. 1214–1217, 2011.

V. P. Nguyen, T. Rabczuk, S. Bordas, and M. Duflot, “Meshless methods: a review and computer implementation aspects,” Mathematics and Computers in Simulation, vol. 79, no. 3, pp. 763–813, 2008.

R. Han, J. Liang, X. Qu et al., “A source reconstruction algorithm based on adaptive hp-FEM for bioluminescence tomography,” Optics Express, vol. 17, no. 17, pp. 14481–14494, 2009.