A Feasibility Study: Operator Splitting for Solving Anisotropic Problem

The electroencephalography (EEG) is a non-invasive technique to study electrical brain activity (while brain is performing a cognitive task). The electrical brain activity is a complex process of electrical propagation because the brain structure is an incredibly complex structure. This complex structure leads to different conductivity property in term of its magnitude and orientation, called anisotropic conductivity. Using Maxwell’s equations, the electrical brain activity has been studied intensively. For simpliﬁcation, the quasistatic Maxwell’s equations are used to model the electrical brain activity and it leads to deal with a Poisson’s equation. In this research, a feasibility study of using a new method, called Operator Splitting Method (OSM), to solve anisotropic 2-Dimensional (2D) Poisson’s equation is performed. A freeware of ﬁnite element method (FEM), FreeFEM++, is employed to build matrices used in the OSM algorithm. The OSM algorithm which is written in Matlab is then tested to solve anisotropic 2D Laplace’s equation and anisotropic Poisson’s equation with dipolar source. Afterwards, the OSM solutions are validated by using exact solution and direct numerical solution. By using L2-Error Norm, the convergence rate of the OSM algorithm is then analyzed. Some numerical experiments have been performed to test the performance of the OSM algorithm. The OSM solution of anisotropic 2D Laplace’s equation coincide with the exact and direct numerical solution of the problem. For anisotropic 2D Poisson’s equation with dipolar source, some similar results has been obtained too. The pattern of the OSM solutions are similar to the pattern of direct numerical solutions of the problem. The results arise a hope to attempt implementing the OSM algorithm for more complex problem such as a realistic human head model.


Introduction
The electroencephalography (EEG) is a non-invasive technique to study electrical brain activity (while brain is performing a cognitive task). The electrical brain activity is a complex process of electrical propagation because the brain structure is an incredibly complex structure. This complex structure leads to different conductivity property in term of its magnitude and orientation, called anisotropic conductivity. Using Maxwell's equations, the electrical brain activity has been studied intensively. In order to implement the Maxwell's equations for EEG modelling, the time-derivatives the Maxwell's equations could be omitted, that is called quasistatic approximation, and it leads to deal with a Poisson's equation as shown by Equation 1. This approximation is valid for EEG modelling because the frequency of EEG on the human head is very weak, lower than 10 2 Hz [1].
∇.(σ∇V ) = ∇.J p (1) where V is the electric potential, J p is the primary cur-rent, and σ is the conductivity properties. In EEG human head modelling, the conductivity properties of equation 1 must be specified. This property is related to the structure of human head. The structure of human head is an incredibly complex structure because it is a heterogeneous structure with different tissue layers and each tissue layer has particular property as shown by Figure 1 [2] [3]. Thus, it is a challenge for scientist to deal with this issue. For a realistic human head model, it is very difficult to solve the model which consider the complexity of human head structure [4] [5] [6]. Generally, this model is called anisotropic human head model. The idea of this report is to propose a new method, called Operator Splitting Method (OSM), to solve anisotropic problem. The method proposed is expected to be able to solve anisotropic problem easily. In this report, a feasibility study of using OSM to solve a simple anisotropic problem is performed. Hopefully, this feasibility study can motivate researchers to implement OSM in a realistic human head model. Many researchers have tried to develop such a human head model for solving forward problem of EEG source localization. The forward problem of EEG aims at computing the electrical field produced by a known primary current, in a known geometry. The boundary ele-ment method (BEM) and the finite element method (FEM) can be implemented for solving the problem. Both of these numerical approximation methods, BEM and FEM, have either advantages or disadvantages. Using BEM will allow to carry out simulations at low computational cost because by this approach, the domain decomposition or meshing is only applied on the boundary as shown in Figure 2a.
The disadvantage of BEM is in non-realistic assumptions used (all tissues are homogenous and isotropic between boundaries) [7]. While using FEM requires high computational cost to perform simulation because the whole computation domain must be decomposed or meshed as shown by Figure 2b. It means that the number of unknown vertices in using FEM is much higher than using BEM. However, FEM approach is much more realistic than BEM approach because anisotropic material parameters could be incorporated in using FEM. In this research, using OSM, a simple anisotropic problem is tried to be solved by considering its isotropic model and introducing a correction term, called lineic. Therefore, the solution of OSM is compared to the exact and direct numerical solution of this anisotropic problem. In 2D domain of xy -coordinates , the conductivity operator in equation 1 and 2, σ , can be presented as a matrix:

Experimental / Theoretical Method
where σ x is the conductivity in x -direction ,and σ y is the conductivity in y -direction. If σ x = σ y , then σ is an isotropic conductivity, σ iso . When σ x φ σ y φ 0 or 0 π σ x π σ y , then σ is an anisotropic conductivity, σ anis . Thus, in xy -coordinates, the 2D Poisson's and Laplace's equation can be stated by equation 4 and 5, respectively.
In order to implement OSM in the problem of interest, the conductivity operator in anisotropic 2D Poisson's and anisotropic 2D Laplace's equation is treated as the summation of an isotropic conductivity, σ iso , and a lineic or correction conductivity, σ lin , as follows: In the xy -coordinates, if the assumption of σ x φ σ y is taken, equation 6 can be represented as a matrix addition as follows: Considering equation 1 and 2 as 2D anisotropic problem in xy -coordinates and substituting equation 6 to both equations, if σ y = 1 , equation 1 and 2 can be re-arranged as equation 8 and 9, respectively.
The value of σ x−lin depends on how stronger σ x , the conductivity in x -direction, to σ y , the conductivity in ydirection. Solving Poisson's equation using FEM will yield a linear system AV = B. Thus, if equation 8 and 9 are solved using FEM, the linear system will be produced as well. Ba-sically, the linear system AV = B is built as follows: matrix A is built from the left hand-side of the equations, while matrix B is built from the right hand-side of the equations. It can be seen that using OSM, there are two terms in the left hand-side of the equations. It means that there are two matrices A that will be produced. The first matrix A is built from , both terms will be discretized and the elements of the matrices will be generated. Afterwards, the OSM algorithm will be developed by using those matrices produced. In this research, the OSM algorithm proposed is as follows: where λ is the relaxation coefficient; constant and n is the iteration number. Introducing a relaxation coefficient, λ, in the algorithm proposed is designed to have a hope of speeding-up the convergence of the value of iteration. The OSM algorithm which is written in Matlab is then tested to solve an anisotropic Laplace's and Poisson's equation with dipolar source. Afterwards, the OSM solutions are validated by using exact solution and direct numerical solution. By using L2-Error Norm, the convergence rate of the OSM algorithm is then analyzed. Let us recall the original problem, Poisson's equation (equation 1), where σ x is sigma times of σ y as follows: The boundary condition of equation 13 will be modified as that Neumann boundary value will be equal to zero, due to the hypothesis that states no current flows outward of the human head is used. The problem is solved in 2D rectangular domain with homogeneous Neumann condition on all boundaries, as shown in Figure 3. The Operator Splitting solution will be compared to the direct numerical solution by using the L2-Error Norm. The L2-Error Norm in solving anisotropic 2D Laplace's and Poisson's equation with dipolar source is formulated as: where V direct is the direct numerical solution of the problem.

The Numerical Solution of Anisotropic 2D Laplace's Equation
Before dealing with anisotropic 2D Poisson's equation, the work will be started by solving anisotropic 2D Laplace's equation numerically and then validate its numerical solution with the analytical solution. First of all, the isotropic 2D Laplace's equation will be solved by using FEM. Let us recall 2D Laplace's equation :  is an anisotropic problem in which σ x = 3σ y . On the boundaries, there are two types of the boundary condition (BC) that are used. They are Dirichlet BC in which the solution is fixed on the boundary and Neumann BC where its normal derivative is fixed on the boundary. In using FEM to solve the problem, a free software FreeFEM++ is employed to mesh the 2D rectangular domain and build the linear system of the problem, AV = B. It means that FreeFEM++ will build both matrices, matrix A and B. Furthermore, a Matlab script will be used to solve this linear system and display the numerical solution in graphical representation.
The dimensions of 2D rectangular domain are 0.1 × 0.7, decomposed using a triangle meshing into 1852 nodes or vertices and 3462 triangles as shown by figure 5. The basis function is choosed to build the linear system. The direct numerical solution of anisotropic Laplace's quation in 2D rectangular domain is shown by figure 6.
In using Operator Splitting Method to solve the problem, a free software FreeFEM++ is employed to build matrix A iso , A lin , and B which will be implemented in the OSM algorithm. The OSM gives a solution of the problem as displayed by figure 7. This solution is obtained by parameters setting of λ, relaxation coefficient = 0.1, n, number of iteration = 30. In order to validate the numerical solution of anisotropic 2D Laplace's equation, the analytical solution of anisotropic 2D Laplace's equation will be obtained. For obtaining the analytical solution, the Separation of Variables method is employed. Figure 8 shows the analytical solution of anisotropic 2D Laplace's equation.   Figure 9 and 10 show the comparison between exact and direct numerical solution and the comparison between exact and Operator Splitting solution anisotropic 2D Laplace's equation, respectively. It can be seen that both direct numerical and Operator Splitting solution graphics coincide with the exact solution. It means that either direct numerical scheme and Operator Splitting algorithm used well suit to solve anisotropic 2D Laplace's equation. Thus, the direct numerical scheme and Operator Splitting algorithm can be proposed to solve the next problem, anisotropic 2D Poisson's equation with dipolar source. Particulary for the Operator Splitting algorithm proposed, the convergence rate of this algorithm can be exploited in order to know how well this algorithm works. For this purpose, the L2-Convergence Norm will be used to analyze the convergence for each n -th iteration. Figure 11 shows the convergence rate of the Operator Splitting algorithm proposed. In general, the algorithm proposed works well because the convergence trend is decrease along the iteration.  In evaluating L2-Error Norm, the relaxation coefficient, λ, is varied in order to know the best setting of this parameter in the OSM algorithm. For the number of itera-tion of 30, the L2-Error Norm is shown by figure 12 and 13. While table 1 gives the information of the minimum L2-Error Norm obtained for each λ.  Figure 12 shows that from λ = 0.03 to λ = 0.4, the higher lamda used, the faster the L2-Error Norm obtained. This fact is also showed in table 4.1. While for λ = 0.5 and λ = 0.6, they have similar trend with λ = 0.2 and λ = 0.07, respectively. In figure 13, the trend line of L2-Error Norm obtained by λ = 0.8 is basically similar to the trend line of using λ = 0.7 and λ = 0.9. It can be seen that using λ higher than 0.6, the L2-Error Norm will not converge. Thus, it can be concluded the maximum relaxation coefficient allowed in the OSM algorithm is about 0.6.

The Numerical Solution of Anisotropic 2D Poisson's Equation
In this part, the anisotropic 2D Poisson's equation with dipolar source will be solved by using Operator Splitting Method. For an experiment, the anisotropic 2D Poisson's equation with dipolar source is specified where σ x = sigma · σ y in which sigma = 8. The solution is obtained by parameters setting of λ, relaxation coefficient = 0.1, n, number of iteration = 30, an initial solution for all vertices = 0 and shown by figure 14. By modifying sigma from 2 to 10, the computation parameters for obtaining a proper solution are shown by table 2.
In solving anisotropic 2D Poisson's equation with dipolar source, the Operator Splitting solution will be validated by using the direct numerical solution. The direct numerical scheme used, FEM, is very well-known technique and it has been used in wide area. Hence, the Operator Splitting solution will be compared to the direct numerical solution by using the L2-Error Norm. Figure  15 shows the comparison between the OSM and direct numerical solution for sigma = 0.8.

Conclusions
The numerical experiments show that the Operator Splitting Method proposed is able to solve either anisotropic 2D Laplace's equation or anisotropic 2D Poisson's equation with dipolar source. In solving anisotropic 2D Laplace's equation, the OSM solution coincides with the exact and direct numerical solution of the problem. While for solving anisotropic 2D Poisson's equation with dipolar source, the pattern of the OSM solutions are similar to the pattern of the direct numerical solutions. Indeed, some of OSM solutions almost coincide with the direct numerical solution.
In solving anisotropic 2D Poisson's equation with dipolar source, the maximum number of iteration is specified equal to 30. This parameter setting is based on some facts that a setting of number of iteration more than 30 will not lead a convergence iteration. Furthermore, for determining computation parameters, the relaxation coefficient and number of iteration, to obtain "the best" solution of the problem, the minimum L2-Error Norm obtained and its number of iteration are consider. Sometime, there is a trade-off between those values in choosing a proper computation parameters because of a small differences of those values.
The performance of the OSM algorithm proposed works well in solving either anisotropic 2D Laplace's equation or anisotropic 2D Poisson's equation with dipolar source. It arises a hope or motivation to attempt implementing the OSM algorithm for more complex problem such as a realistic human head model. Applying the OSM algorithm in a such complex domain, it seems that the most difficult task is building matrix that will be implemented to the algorithm proposed. Thus, in the future, this challenge is a very interesting research topic for researchers. And hopefully, a new method, Operator Splitting Method, will be really able to solve a real problem of EEG modelling.