# FEM method for the EEG forward problem and improvement based on modification of the saint venant's method

Takfarinas Medani, David Lautru, Denis Schwartz, Zhuoxiang Ren, Gérard Sou

<figure><img src="image_1.png" /></figure>

## To cite this version:

Takfarinas Medani, David Lautru, Denis Schwartz, Zhuoxiang Ren, Gérard Sou. FEM method for the EEG forward problem and improvement based on modification of the saint venant's method. Progress In Electromagnetics Research, 2015, 153, pp.11-22. 10.2528/PIER15050102. hal-01222242

HAL Id: hal-01222242

https://hal.sorbonne-universite.fr/hal-01222242

Submitted on 18 Jan 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# FEM Method for the EEG Forward Problem and Improvement
Based on Modification of the Saint Venant's Method

Takfarinas Medani¹, *, David Lautru², Denis Schwartz³,
Zhuoxiang Ren¹, and Gerard Sou¹

**Abstract**—The finite-element method (FEM) is applied to solve the EEG forward problem. Two issues related to the implementation of this method are investigated. The first is the singularity due to the punctual dipole sources and the second is the numerical errors observed near the interface of different tissues. To deal with the singularity of the punctual dipole sources, three source modeling methods, namely, the direct, the subtraction and the Saint Venant's methods, are examined. To solve the problem of numerical instability near the interface of different tissues, a modification on the Saint Venant's method is introduced. The numerical results are compared with analytical solution in the case of the multilayer spherical head models. The advantages of the proposed method are highlighted.

## 1. INTRODUCTION

Understanding the brain and its functions in vivo is a predominant research field. To this end, many devices are developed, and among them the electroencephalography (EEG) is a non-invasive functional imaging technique. The EEG directly measures the electric brain activity on the scalp with high temporal resolution. Its main purpose is to record the brain activity during a specific task and then to localize the sources of the recorded signals (neurons populations) inside the brain. This tool is used both in clinical practice and cognitive research areas in neuroscience [1, 2]. The macroscopic neuronal activity is assumed to be modeled by electrical current dipole sources [2]. The process of source localization involves the resolution of an inverse problem which requires several models to solve the associated forward problem. The forward problem allows the prediction of the potential on the scalp generated from a given source inside the brain. Its resolution leads to constructing the so called “lead field matrix” [2, 3], an operator describing the relation between the source parameters, the head parameters and the predicted signals on the EEG electrodes. The accuracy of the source localization depends largely on the accuracy of the forward problem’s solution, which is related to the choice of the head’s model, source’s model and computational method.

Currently, in practice, for the sake of computational simplicity and rapidity, source's localization software mostly uses the multilayer spherical head model with an analytical solution, or the boundary element method (BEM) with a numerical solution to solve the forward problem [1, 3, 4, 16, 21], where the brain tissues are treated as homogenous and isotropic. The improvement of EEG imaging accuracy requires more realistic head model that takes into account the whole head with more realistic properties [5, 6, 8]. For these reasons, methods based on a volume discretization, such as the finite element method (FEM) or the finite difference method (FDM), have attracted more attention in recent

Received 1 May 2015, Accepted 22 June 2015, Scheduled 12 October 2015

* Corresponding author: Takfarinas Medani (takfarinas.medani@upmc.fr).

¹ The Electronics and Electromagnetism Laboratory, UPMC Univ. Paris 06, UR2, L2E F-75005, Paris, France. ² The Energy Laboratory Energetics Mechanical Electromagnetism, LEME EA 4416 Univ. Paris Ouest Nanterre La Defense, Ville d'Avray, France.

³ The Brain & Spine Institute (ICM), UMR 7225/U 1127, UPMC/CNRS/INSERM, MEG-EEG-CENIR, Pitié-Salpêtrière Hospital, Paris, France.
years [7–9, 13]. However, with a FEM discretization, the punctual source introduces a numerical singularity which impacts on both the accuracy and the performance of the FEM forward solution. The direct method, subtraction method and saint Venant’s method [4, 5, 7, 9, 10, 14, 15, 20] are three approaches which can be used to improve the behavior of the FEM in this context. In the literatures, several studies investigate and compare the accuracy of these methods [9, 10, 14–17]. However, the previous studies observed remaining numerical instability in the case of sources located near the interfaces of different layers. In this study, these three aforementioned methods are investigated and applied in the case of multilayer spherical head models where the current dipole sources are placed in different positions and orientations inside the brain’s layer. To deal with the problem of the numerical instability near the interfaces of different tissues, a modification on the Venant’s method is introduced and evaluated. The resulting surface potential is compared with those obtained from an available analytic solution [11].

## 2. THE EEG FORWARD PROBLEM

The EEG forward problem is the computation of the electric potential V on the scalp from a given source inside the brain. This calculation leads to solving the Poisson’s equation, which needs models of the source (brain activities) and the head.

### 2.1. Source Model

The primary source of the EEG signals is a large number of active pyramidal neurons, consisting of around $10^5$ to $10^7$ cells, laying in a small volume (called macro-column) and located in the cortex, i.e., the grey matter. The pyramidal neurons are aligned perpendicularly to the cortex’s surface. They are synchronously activated during a small period of time (about 10 ms) then their individual electrical contributions can be summed and reach detectable levels. This activity can be modelled by a current dipole [2]. In bioelectricity [2, 7] the dipole can be represented by a pair of monopoles having equal magnitude I but opposite polarity, separated by an infinitesimal distance l as expressed in the relation (1) and shown in Fig. 1.

$$ \vec{j}_p = I [\delta(\vec{r} - \vec{r}_{\text{source}}) - \delta(\vec{r} - \vec{r}_{\text{sink}})] \quad (1) $$

where $\vec{j}_p$ is the primary current density composed of the monopole source at $\vec{r}_{\text{source}}$ position and the monopole sink at $\vec{r}_{\text{sink}}$ position, and $\delta$ is the Dirac distribution. In electromagnetic, this relation can be represented by the dipolar moment vector $\vec{M}$, oriented from the sink to the source and applied at the middle of the segment $[\vec{r}_{\text{sink}}, \vec{r}_{\text{source}}]$ at $\vec{r}_0$:

$$ \vec{j}_p = I \vec{l} \delta(\vec{r} - \vec{r}_0) = \vec{M} \delta(\vec{r} - \vec{r}_0) \quad (2) $$

The associated current density in the conductor $\vec{j}_s$, also called the volume or secondary current is given by the Ohm’s law:

$$ \vec{j}_s = \sigma \vec{E} \quad (3) $$

where $\vec{E}$ denotes the electric field in the volume conductor, and $\sigma$ is the electric conductivity. The total current $\vec{j}$ in the volume is:

$$ \vec{j} = \vec{j}_p + \vec{j}_s \quad (4) $$

The bioelectric field’s behaviour in the head can be expressed by Maxwell’s equations. The frequencies of EEG signals are up to 100 Hz, hence, the propagation and time variation can be ignored [2]. Therefore the quasi-static approximation is used:

$$ \operatorname{div} \vec{j} = 0 \quad (5) $$

$$ \vec{E} = -\operatorname{grad} V \quad (6) $$

These equations lead to the Poisson’s equation:

$$ \operatorname{div} [\sigma \operatorname{grad} V] = \operatorname{div} \vec{j}_p = J_p \quad (7) $$
<figure><img src="image_3.png" /><figcaption>Figure 1. Dipole moment's representation.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. The multi-layer concentric sphere.</figcaption></figure>

## 2.2. Head model

EEG source analysis usually uses a spherical model to approximate the human head. It consists of N layers, representing the different tissues with the radii $r_i$ and the conductivity $\sigma_i$, $i = 1, ..., N$ for the ith layer, an example is shown in the Fig. 2. The solution of electric potential is evaluated on the outermost surface $S_N$. This model is used to study the forward and inverse EEG problem in several studies [1-4]. Although, the spherical model is not accurate because it does not represent the real human head, however, it plays an important role in the research since it can provide analytical solution [11] and serve as reference to validate the numerical methods.

In this paper two spherical models of the head are used. The first is built with four concentric layers, noted as $hm1$, and the second is represented in three layers and noted as $hm2$. Starting from the inside, for the model $hm1$, the layers represent, respectively, the brain, the cerebro-spinal fluid (csf), the skull and the scalp. The outside radius of each layer are [63, 65, 71, 75] in mm and the conductivities are [0.33, 1.0, 0.0042, 0.33] in $Sm^{-1}$, respectively. For the second model $hm2$, the layers are the same as previous one but without the csf's compartment, the outer radius of each layer are [63, 71, 75] (mm), respectively.

## 3. FEM FORMULATION OF THE FORWARD PROBLEM

To apply the FEM, the entire volume $\Omega$ representing the head is divided into tetrahedral elements; the nodal potential $V$ is expressed as:

$$V(x, y, z) = \sum_{i=1}^{n} V_i \varphi_i(x, y, z) \qquad (8)$$

where $\varphi_i(x, y, z)$ denotes a set of test functions also called basis functions, $V_i$ is the scalar potential associated to the node $i$, $n$ is the number of vertices (FEM nodes). On the outer surface, the following Neumann boundary condition is imposed:

$$\sigma \overrightarrow{grad} V \cdot \vec{n} = 0 \qquad (9)$$

which means that the normal component of the current density vanishes, $\vec{n}$ is the normal vector to the surface $S_N$.

To derive to the finite element formulation, the Galerkin approach is applied and the 'weak formulation' is obtained:

$$\int_{\Omega} \overrightarrow{grad} \varphi \cdot \sigma \overrightarrow{grad} V d\Omega = \int_{\Omega} \varphi \operatorname{div} \vec{j}_p d\Omega \qquad (10)$$

Substituting (8) in (10) leads to the following linear matrix system:

$$KV_i = J \qquad (11)$$

where $J$ is the source term obtained by integration of the right hand side (rhs) of Equation (10), and $K$ is the stiffness matrix.

In the implementation of the FE approach on the EEG forward problem, there are two issues need to be investigated. The first is the treatment of the singularity caused by the puncture source. The second
is the important numerical instabilities observed for the sources situated near the interface of different conductivities [5, 9, 10]. In the following sections, three dipole source modelling methods, namely, the direct method, subtraction method and Saint Venant's method, are respectively, described. The Saint Venant's method is further investigated to deal with the numerical problems near the interfaces, where a modification on the monopoles' distribution, based on the Saint Venant Condition [14, 17, 20], is introduced and tested.

# 4. TREATMENT OF THE SOURCE'S SINGULARITY

The numerical singularity of the point source affects the accuracy of the final results. In this section the direct method, the subtraction method, and the Saint Venant's method dealing with the singularity are presented and evaluated.

## 4.1. Direct Method

This method is well described in [4, 9, 10, 14]. It is formulated by integrating the dipole source with the element basis functions $\varphi_i$. Its formulation is expressed in the following:

$$ \int_{\Omega} \overrightarrow{\text{grad}} \varphi_i \cdot \sigma \overrightarrow{\text{grad}} V d\Omega = - \int_{\Omega} \vec{j}_p \overrightarrow{\text{grad}} \varphi_i d\Omega \quad (12) $$

The right hand side (rhs) of the Poisson's equation is defined only on nodes of the element containing the dipole and set to zero otherwise as shown in Fig. 3. As can be easily observed, the rhs of (12) is highly impacted by the shape of the element where the source is located.

## 4.2. Subtraction Method

The idea is to split the total potential $V$ into a singularity potential $V^\infty$ and a correction potential $V^{\text{corr}}$ [4, 5, 9, 10].

$$ V = V^{\infty} + V^{\text{corr}} \quad (13) $$

The potential $V^\infty$ is analytically calculated for a dipole source at position $\vec{r}_0$ in a homogeneous volume conductor $\Omega^\infty$ with constant isotropic conductivity $\sigma^\infty$. From the Equations (7), (9) and (13), these following relations are obtained:

$$ \text{div} \left[ \sigma \overrightarrow{\text{grad}} V^{\text{corr}} \right] = \text{div} \left[ (\sigma - \sigma^{\infty}) \overrightarrow{\text{grad}} V^{\infty} \right] \quad (14) $$

$$ \sigma \overrightarrow{\text{grad}} V^{\text{corr}} \cdot \vec{n} = - \sigma \overrightarrow{\text{grad}} V^{\infty} \cdot \vec{n} \quad (15) $$

$V^{\text{corr}}$ is numerically calculated by solving (14). Finally the total potential $V$ can be calculated using the relation (13). The right hand side of the Poisson's equation in this method is defined on all nodes of elements that have different conductivity than $\sigma^\infty$.

## 4.3. Saint Venant's Method

This approach is based on the Saint Venant's principle [5, 10, 12], which is originally applied in the elasticity's theory. It proposes to substitute the dipole by a set of monopoles, placed around the dipole's position $\vec{r}_0$. The values of these monopoles are determined to approximate as well as possible the dipole moment $\vec{M}$ by the multi-pole moment $\vec{T}$ given by this relation:

$$ \vec{T} = \int_{\Omega} (\vec{r}_v - \vec{r}_0) J_p d\Omega \quad (16) $$

where $\vec{r}_v$ denotes the position of the monopoles, and $J_p = \text{div} \vec{j}_p$ is the volume current density (Am$^{-3}$). In the case of the FEM, the relation (16) is discretized; and each monopole is placed on a neighbour FE node. Hence the relation (16) yields to:

$$ \vec{T} = \sum_{v=1}^{k} (\vec{r}_v - \vec{r}_0) \vec{j}_{ven} = \sum_{v=1}^{k} \vec{r}_v \delta_{ven} = \bar{R} \vec{j}_{ven} \quad (17) $$
where $k$ is the number of monopoles, $j_{ven}$ the discrete value of the monopole, $\vec{r}_{v0}$ the vector from the node $v$ at position $\vec{r}_v$ to the real dipole position $\vec{r}_0$, and $\bar{R}$ the matrix of distances $\vec{r}_{v0}$. The values $j_{ven}$ are determined by the minimisation of the following expression (18) using the Tikhinove regularization:

$$ \tau = \arg \min_{j_{ven}} \left[ \left\| \vec{M} - \bar{R} j_{ven} \right\|_2^2 + \lambda \left\| \bar{W} j_{ven} \right\|_2^2 \right] \quad (18) $$

with $\|u\|_2$ the discrete $L^2$ norm. The additional of the second term in (18) constrains the solution, depending on the choice of the regularization matrix $\bar{W}$ and the corresponding parameter $\lambda$ [5, 12, 14].

Solving the above minimization system with the first differentiation with respect to $j_{ven}$ gives the following expression:

$$ j_{ven} = \left[ \left[ \sum_{d=x,y,z} (\bar{R})_d^T \bar{R}_d + \lambda (\bar{W})_d^T \bar{W}_d \right]^{-1} \left[ \sum_{d=x,y,z} [(\bar{R})_d^T M_d] \right] \right] \quad (19) $$

The indexes $d$ and $T$ refer respectively to the direction $(x, y$ or $z)$ and the transposed matrix.

The divergence of the primary current $\vec{j}_p$ using the discrete load distribution can be expressed by the following:

$$ \mathrm{div} \vec{j}_p = \sum_{v=1}^{k} j_{ven} \delta(\vec{r} - \vec{r}_v) \quad (20) $$

$$ \int_{\Omega} \varphi \mathrm{div} \vec{j}_p d\Omega = \sum_{v=1}^{k} j_{ven} \int_{\Omega} \varphi \delta(\vec{r} - \vec{r}_v) d\Omega \quad (21) $$

here the term $\delta(\vec{r} - \vec{r}_v)$ is equal to 1 when the position of the node $\vec{r}$ corresponds to neighboring FE node $\vec{r}_v$ and to zero otherwise, and the function $\varphi$ is either one or null. The rhs for this method, in our study, is defined only on the neighbouring FE nodes around the dipole's position as shown in Fig. 4.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The arrow is the dipole; the dots are the nodes of the tetrahedral element, where the source term is distributed with the direct method.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Example on 2D, the arrow is the dipole; the dots are the nodes where the monopoles are distributed with the Venant's method.</figcaption></figure>

## 5. RESULTS IN SPHERICAL HEAD MODEL

The head model is discretized by a regular tetrahedral mesh of 115734 elements and 20320 nodes using iso2mesh [18]. We compute the Joe-Liu mesh parameter [19] which measures the mesh quality, the average value of quality is 84%.

An equivalent current dipole has been selected to represent the electrical activity in the brain region, with the moment magnitude of $M = 1$ nA·m placed on the z-axis, with various positions and orientations.

The potential distribution is evaluated on the 1425 surface's nodes. These nodes are localized by their elevation angle $\theta \in [-\frac{\pi}{2}, +\frac{\pi}{2}]$ and the azimuthally orientation $\vartheta \in [0, 2\pi]$. The $\theta = 0$ corresponds
<figure><img src="image_7.png" /><figcaption>**Figure 5.** View of the mesh; the $\theta$ and $\emptyset$, (a) a radial dipole and (b) a tangential dipole.</figcaption></figure>

to the plane defined by the *x* and *y* axes and $\emptyset = 0$ to the plane defined by the *x* and *z* axes. Fig. 5 illustrates a view of the mesh, the representation of radial and tangential dipoles and an example of $\theta$ and $\emptyset$.

**5.1. Potential Distribution on the Outer Surface of the Spherical Model**

Figure 6 shows the potential distribution on the nodes of the outer sphere, obtained with the Saint Venant's method, and the numerical (*Vfem) and analytical (°Vana) solutions are plotted on the same graph as shown.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Analytical (°Vana) and numerical solution (*Vfem) on all surface's nodes versus the theta angle (radian), (b) shows the result for a centered radial dipole, (a) those for a centered tangential dipole.</figcaption></figure>

For the radial and tangential dipole, the numerical solution is identical to the analytical solution on each node; a similar potential's distribution is obtained with the others methods (direct and subtraction).

**5.2. Comparison and Discussion of the Results**

The errors between the numerical and analytical solutions are evaluated with different measures of similarity [2, 5–7]: the Relative Difference Measure (RDM) and Magnitude factor (MAG) are used,
respectively, for changes in the topography and the magnitude between the two solutions. The correlation (COR) is used as a measure of similarity between two data sets.

$$ \text{RDM}(V^{ana}, V^{fem}) = \left\| \frac{V_i^{fem}}{\|V_i^{fem}\|} - \frac{V_i^{ana}}{\|V_i^{ana}\|} \right\|_2 \qquad (22) $$

$$ \text{MAG}(V^{ana}, V^{fem}) = \frac{\|V_i^{fem}\|_2}{\|V_i^{ana}\|_2} \qquad (23) $$

$$ \text{COR}(V^{ana}, V^{fem}) = \frac{\sum_{i=1}^{N} (V_i^{fem} - V_{mean}^{fem}) (V_i^{ana} - V_{mean}^{ana})}{\sqrt{\sum_{i=1}^{N} (V_i^{fem} - V_{mean}^{fem})^2} \sqrt{\sum_{i=1}^{N} (V_i^{ana} - V_{mean}^{ana})^2}} \qquad (24) $$

where $V_i^{fem}$ and $V_i^{ana}$ indicate the numerical and the analytical solutions respectively at the node i. The notations *mean* indicates the mean value of the associated solution. These criterions are evaluated for various dipole positions and orientations. Table 1 (and Table 2) present the results obtained for a radial dipole (and a tangential dipole) with various eccentricities along the z-axis.

For the first head model *hm1*, the interface brain-csf is located at *r* = 63 mm. According to the above results, we note that for the source located deep inside the brain layer, at *r* < 55 mm, the subtraction method and the Saint Venant's method show the best RDM, both for the tangential and the radial orientations. For these positions, the direct approach has the largest errors. The *MAG*

<table><caption>**Table 1.** Results for radial dipole on different positions along the z-axis.</caption><thead><tr><th rowspan="2">hm1</th><th rowspan="2">Method</th><th colspan="6">Position</th></tr><tr><th>0 mm</th><th>30 mm</th><th>55 mm</th><th>62 mm</th><th>62.5 mm</th><th>63 mm</th></tr></thead><tbody><tr><td rowspan="3">RDM (%)</td><td>Direct</td><td>2.05</td><td>1.23</td><td>2.48</td><td>1.64</td><td><strong>1.34</strong></td><td><strong>6.29</strong></td></tr><tr><td>Subtraction</td><td>0.02</td><td>0.18</td><td>1.86</td><td><strong>16.15</strong></td><td><strong>16.80</strong></td><td><strong>17.11</strong></td></tr><tr><td>St Venant</td><td>0.24</td><td>0.53</td><td>0.45</td><td>0.72</td><td><strong>7.02</strong></td><td><strong>6.79</strong></td></tr><tr><td rowspan="3">MAG</td><td>Direct</td><td>0.99</td><td>1.00</td><td>0.99</td><td>1.01</td><td><strong>1.00</strong></td><td><strong>0.74</strong></td></tr><tr><td>Subtraction</td><td>1.00</td><td>1.00</td><td>1.07</td><td><strong>24.55</strong></td><td><strong>31.18</strong></td><td><strong>29.28</strong></td></tr><tr><td>St Venant</td><td>1.00</td><td>1.00</td><td>1.01</td><td>0.79</td><td><strong>10.96</strong></td><td><strong>13.31</strong></td></tr><tr><td rowspan="3">COR (%)</td><td>Direct</td><td>99.88</td><td>99.95</td><td>99.84</td><td>99.93</td><td><strong>99.95</strong></td><td><strong>98.05</strong></td></tr><tr><td>Subtraction</td><td>100</td><td>99.99</td><td>99.90</td><td><strong>93.49</strong></td><td><strong>92.98</strong></td><td><strong>92.75</strong></td></tr><tr><td>St Venant</td><td>99.99</td><td>99.99</td><td>99.99</td><td><strong>99.98</strong></td><td><strong>98.83</strong></td><td><strong>98.91</strong></td></tr></tbody></table>

<table><caption>**Table 2.** Results for tangential dipole on different positions along the z-axis.</caption><thead><tr><th rowspan="2">hm1</th><th rowspan="2">Method</th><th colspan="6">Position</th></tr><tr><th>0 mm</th><th>30 mm</th><th>55 mm</th><th>62 mm</th><th>62.5 mm</th><th>63 mm</th></tr></thead><tbody><tr><td rowspan="3">RDM (%)</td><td>Direct</td><td>0.36</td><td>1.71</td><td>0.99</td><td>1.19</td><td><strong>1.43</strong></td><td><strong>4.12</strong></td></tr><tr><td>Subtraction</td><td>0.02</td><td>0.15</td><td>1.13</td><td><strong>8.62</strong></td><td><strong>4.57</strong></td><td><strong>3.52</strong></td></tr><tr><td>St Venant</td><td>0.18</td><td>0.56</td><td>1.10</td><td>0.76</td><td><strong>34.10</strong></td><td><strong>22.16</strong></td></tr><tr><td rowspan="3">MAG</td><td>Direct</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td><td><strong>0.99</strong></td><td><strong>1.99</strong></td></tr><tr><td>Subtraction</td><td>1.00</td><td>1.00</td><td>1.01</td><td><strong>1.22</strong></td><td><strong>1.60</strong></td><td><strong>2.17</strong></td></tr><tr><td>St Venant</td><td>1.00</td><td>1.00</td><td>1.01</td><td>1.00</td><td><strong>1.53</strong></td><td><strong>1.70</strong></td></tr><tr><td rowspan="3">COR (%)</td><td>Direct</td><td>99.99</td><td>99.94</td><td>99.98</td><td>99.97</td><td><strong>99.95</strong></td><td><strong>98.15</strong></td></tr><tr><td>Subtraction</td><td>100</td><td>99.99</td><td>99.97</td><td><strong>98.95</strong></td><td><strong>98.67</strong></td><td><strong>98.78</strong></td></tr><tr><td>St Venant</td><td>99.99</td><td>99.99</td><td>99.97</td><td>99.98</td><td><strong>81.58</strong></td><td><strong>90.06</strong></td></tr></tbody></table>
error is bounded between 0.99 and 1.01. The variations are not important, and no big differences are noted between the three approaches. The coefficient correlation (COR) shows good values for the three methods for these source locations.

When the source is located between 55 mm ≤ *r* ≤ 62 mm, the subtraction method loses its accuracy, the RDM rises, whereas the Venant's method and the direct methods remain stable. The *MAG* and *COR* follow the same trend as the *RDM*.

For the most eccentric source's positions, located between 62 mm < *r* < 63 mm the Saint Venant's method also loses its performance. Moving the dipole from 62 mm to 62.5 mm induces RDM jumps from 0.72% to 7.02% in the case of a radial orientation. The jump of errors is observed for the tangential orientation as well. A relative stability for the direct method is observed.

As regards to the sources on the interface, *r* = 63 mm, all the three approaches give bad results. Similar observations are highlighted in the previous publications, the numerical errors rise when the dipole approaches the interface [5, 7, 10, 17].

The majority of the sources are located in the cortex, i.e., the real position of the dipole is often close to the interface, the accuracy of results in this situation is very important. However, in these positions important errors are observed for all the studied methods. To deal with this problem, in this article, we introduce an original modification on the saint Venant's method for the dipoles located near the interface, and this modification is based on the "Venant's Condition" described in the following section.

# 6. THE MODIFIED SAINT VENANT'S METHOD

Based on the "Venant Condition", namely no sources should be placed in other than grey matter compartments [14, 16, 17, 20], and due to the stability of its results as shown in Table 1 and Table 2. The monopoles distribution is investigated.

## 6.1. Saint Venant's Monopoles Distribution

In the traditional Saint Venant's approach, the monopoles are distributed on all the neighbouring nodes around and close to the dipole's position, as shown in Fig. 5. For the dipole located near the interface, some of the neighbour monopoles could be located on different layers. From the previous results shown in Table 1 and Table 2, for the dipole's positions at *r* = 62 mm and *r* = 62.5 mm, high jumps on all criterions of errors are observed for the Saint Venant's method. In order to investigate the causes of the sudden increase in error we analysed the locations of the monopoles', and Table 3 shows the result.

<table><caption>**Table 3.** Monopoles' location near the interface.</caption><thead><tr><th>Dipole position</th><th>r = 62 mm</th><th>r = 62.5 mm</th></tr></thead><tbody><tr><td>Monopoles' localisation</td><td></td><td></td></tr><tr><td>In the brain's layer (σ<sub>br</sub> = 0.33 Sm<sup>-1</sup>)</td><td>21</td><td>20</td></tr><tr><td>In the csf's layer (σ<sub>cs</sub> = 1 Sm<sup>-1</sup>)</td><td>9</td><td>9</td></tr><tr><td>In the skull's layer (σ<sub>sk</sub> = 0.0042 Sm<sup>-1</sup>)</td><td>3</td><td>8</td></tr><tr><td>In the scalp's layer (σ<sub>sc</sub> = 0.33 Sm<sup>-1</sup>)</td><td>0</td><td>0</td></tr></tbody></table>

From this table, we observe that there is a set of monopoles located outside of the brain's layer. This distribution does not fulfil the "Venant's Condition" and this could be one of the reasons for the rising on the error near the interface. In addition, moving the dipole from 62 mm to 62.5 mm implies high variation on the monopoles number on the skull's layer; this number passes from 3 monopoles to 8 monopoles. We suppose that the high difference between the conductivity of the brain and the skull associated with the high number of the monopoles in the skull's layer is the cause on the increased errors. Therefore, we propose to calculate the numerical solution by considering only the monopoles inside the brain layer. We remove all the monopoles outside the brain's layer, as shown in Fig. 7.

The mathematical procedure to calculate the monopoles' values is the same as explained in 4.3. The difference is on the number of the monopoles, their locations and values.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) is the original Saint Venant distribution, (b) is the modified distribution including only monopoles inside the brain layer.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Analytical (°Vana) and numerical (*Vfem) potential distribution for a radial dipole near the interface. (a) is the result with all monopoles, (b) those for the monopoles located only in the brain layer.</figcaption></figure>

<table><caption>**Table 4.** Comparison of the results near the interface brain-csf.</caption><thead><tr><th rowspan="3">hm1</th><th>Position</th><th colspan="2">62 mm</th><th colspan="2">62.5 mm</th><th colspan="2">63 mm</th></tr><tr><th>Method</th><th>old</th><th>new</th><th>old</th><th>new</th><th>old</th><th>new</th></tr><tr><th>Orientation</th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td rowspan="2">RDM (%)</td><td>Radial</td><td>0.72</td><td>1.30</td><td><u>7.02</u></td><td>1.19</td><td><u>6.79</u></td><td>3.63</td></tr><tr><td>Tangential</td><td>0.76</td><td>0.91</td><td><u>34.10</u></td><td>1.26</td><td><u>22.16</u></td><td>1.42</td></tr><tr><td rowspan="2">MAG</td><td>Radial</td><td>0.79</td><td>1.00</td><td><u>10.96</u></td><td>1.00</td><td><u>13.31</u></td><td>1.00</td></tr><tr><td>Tangential</td><td>1.00</td><td>0.99</td><td><u>1.53</u></td><td>1.00</td><td><u>1.70</u></td><td>0.99</td></tr><tr><td rowspan="2">COR (%)</td><td>Radial</td><td>99.98</td><td>99.95</td><td><u>98.83</u></td><td>99.92</td><td><u>98.91</u></td><td>99.65</td></tr><tr><td>Tangential</td><td>99.98</td><td>99.98</td><td><u>81.58</u></td><td>99.97</td><td><u>90.06</u></td><td>99.95</td></tr></tbody></table>

## 6.2. Results in Four Concentric Spheres

The modified Venant's method is applied on the *hm1* head model for same dipoles' positions as previous. Fig. 8 shows an example of the potential distribution on the outer surface's node for a radial dipole at the position *r* = 62.5 mm. The abscissa axis represents the *θ* angle which is the same as in Section 5. The first graph shows the result with the standards Venant's method and the second the results with the modified method.

In Table 4 the notations “old” and “new” refer respectively to the standard and the modified Saint Venant’s method.

The results obtained with the modified Saint Venant's method show good improvements on potential distribution for the most eccentric sources. On the other hand, we observe a small rising
on the *RDM* for the sources at *r* = 62 mm but significant improvement on the *MAG*. According to these results, we can say that the error' increase is due to the distribution of the monopoles on several layers with different values of conductivity. Removing the monopoles located outside of the brain layer improves the numerical solution. This modification allows the validation of the “Venant’s Condition” formulated in [14, 17]. To further prove this condition, we apply the same methods on the second head model noted *hm2*.

## 6.3. Results in Three Concentric Spheres Head Model

The head model *hm2* consists of three layers. Compared to *hm1*, the model *hm2* does not have a csf's compartment. The interface is located at *r* = 63 mm, and it is between the brain and the skull. Same simulations as above are evaluated in this model. Table 5 summarizes the errors near the interface and compares the result with the standard and the modified Venant's methods.

<table><caption>**Table 5.** Comparison of the results near the interface brain-skull.</caption><thead><tr><th rowspan="2">hm2</th><th>Position</th><th colspan="2">60 mm</th><th colspan="2">61 mm</th><th colspan="2">62 mm</th></tr><tr><th>Method<br>Orientation</th><th>old</th><th>new</th><th>old</th><th>new</th><th>old</th><th>new</th></tr></thead><tbody><tr><td rowspan="2">RDM (%)</td><td>Radial</td><td>0.65</td><td>0.67</td><td><u>2.39</u></td><td>0.95</td><td><u>1.50</u></td><td>1.65</td></tr><tr><td>Tangential</td><td>0.73</td><td>0.58</td><td><u>38.50</u></td><td>0.67</td><td><u>26.41</u></td><td>1.31</td></tr><tr><td rowspan="2">MAG</td><td>Radial</td><td>1.00</td><td>1.00</td><td><u>10.47</u></td><td>1.00</td><td><u>26.03</u></td><td>1.00</td></tr><tr><td>Tangential</td><td>1.00</td><td>1.00</td><td><u>2.24</u></td><td>0.99</td><td><u>2.14</u></td><td>0.99</td></tr></tbody></table>

<table><caption>**Table 6.** Monopoles' locations versus the dipole position near the interface.</caption><thead><tr><th>Dipole position</th><th>r = 60 mm</th><th>r = 61 mm</th></tr></thead><tbody><tr><th>Monopoles' localisation</th><td></td><td></td></tr><tr><th>In the brain's layer (σ<sub>br</sub> = 0.33 Sm<sup>-1</sup>)</th><td>24</td><td>21</td></tr><tr><th>In the skull's layer (σ<sub>sk</sub> = 0.0042 Sm<sup>-1</sup>)</th><td><u>3</u></td><td><u>12</u></td></tr><tr><th>In the scalp's layer (σ<sub>sc</sub> = 0.33 Sm<sup>-1</sup>)</th><td>0</td><td>0</td></tr></tbody></table>

In Table 5 important variations of the errors between the dipole positions *r* = 60 mm and *r* = 61 mm are observed with the traditional Saint Venant's method, therefore the monopoles' location is investigated in these positions. Table 6 shows the location of the Saint Venant's monopoles for these dipole's positions. The monopoles number on the skull's layer increases from 3 to 12 monopoles, as in the previous case this variation causes the rise in the errors. When these monopoles are removed, the solution is strikingly improved. These results for this model also prove and validate the “Venant’s Condition”.

Furthermore, the simulation' results show that these errors are more important when the difference between the dipole's location's conductivity and the monopole's location's conductivity is higher.

## 6.4. Influence of the Number and the Distribution of Monopoles

In this section, another monopoles' distribution is investigated, this distribution has more monopoles than the distribution used in the previous section. We include not only the neighbours nodes to the dipole, but also those placed close to these neighbours nodes. To make the distinction, we call this new distribution “second cell” noted “2nd” and the original distribution “first cell” noted “1st”. We performed the same simulations as previously described. The following tables summarize the obtained results, respectively, in Table 7 for the *hm1* and in Table 8 for the *hm2*.

With the second cell distribution, the number of the monopoles is increased, the obtained result shows slight improvement compared with the results of the first cell distribution. The results shown in
<table><caption>**Table 7.** Results near/on the interface brain-csf for the two distributions in *hm1*.</caption><thead><tr><th rowspan="2">hm1</th><th>Position</th><th colspan="2">62 mm</th><th colspan="2">62.5 mm</th><th colspan="2">63 mm</th></tr><tr><th>Distribution</th><th>1st</th><th>2nd</th><th>1st</th><th>2nd</th><th>1st</th><th>2nd</th></tr></thead><tbody><tr><td rowspan="2">RDM (%)</td><td>Radial</td><td>1.30</td><td>0.61</td><td>1.19</td><td>0.77</td><td>3.63</td><td>1.04</td></tr><tr><td>Tangential</td><td>0.91</td><td>0.86</td><td>1.26</td><td>1.21</td><td>1.42</td><td>1.12</td></tr><tr><td rowspan="2">MAG</td><td>Radial</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td></tr><tr><td>Tangential</td><td>0.99</td><td>0.99</td><td>1.00</td><td>0.99</td><td>0.99</td><td>0.99</td></tr></tbody></table>

<table><caption>**Table 8.** Results near/on the interface brain-skull for the two distributions in *hm2*.</caption><thead><tr><th rowspan="2">hm2</th><th>Position</th><th colspan="2">62 mm</th><th colspan="2">62.5 mm</th><th colspan="2">63 mm</th></tr><tr><th>Distribution</th><th>1st</th><th>2nd</th><th>1st</th><th>2nd</th><th>1st</th><th>2nd</th></tr></thead><tbody><tr><td rowspan="2">RDM (%)</td><td>Radial</td><td>1.65</td><td>1.11</td><td>2.21</td><td>1.5</td><td>4.39</td><td>1.94</td></tr><tr><td>Tangential</td><td>1.31</td><td>1.24</td><td>1.55</td><td>1.7</td><td>1.87</td><td>1.71</td></tr><tr><td rowspan="2">MAG</td><td>Radial</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td></tr><tr><td>Tangential</td><td>0.99</td><td>0.99</td><td>0.98</td><td>0.98</td><td>0.98</td><td>0.98</td></tr></tbody></table>

Table 7 are better than all the results of Table 4. These outcomes indicate that taking into account more monopoles inside the brain's layer, in the case of high eccentricity, can further enhance the numerical solution with the modified Venant's method.

## 7. CONCLUSION

In this study the finite element method (FEM) is investigated to solve the forward problem and calculate the scalp's electrical potential from the dipole current source. The singularity of the source introduced by the FEM discretization is treated with the direct method, subtraction method and Saint Venant's method. These methods are investigated in multilayer sphere models with an available analytical solution. Results show that all methods give good results in the case of sources located deeply inside the brain's layer. On the contrary, in the case of sources close to the interface of layers with different conductivities, the observed results show numerical instability for all approaches.

To deal with this issue, a modification is introduced to the Saint Venant's method through the implementation of the so called "Venant Condition". This condition states that "no monopole should be placed in other location than grey matter". Therefore, the modification conserves only the monopoles inside the brain's layer in the calculation. The modified Saint Venant's method gives satisfactory results and improves both the topography (with the RDM measure) and the magnitude (the MAG measure) for the sources with highest eccentricities. To further improve the accuracy, a new distribution of the monopoles of the modified Venant's method is introduced. With this addition, better results are obtained compared to the previous methods.

In conclusion, the accuracy of the Saint Venant's method is linked with the monopoles' positions, monopoles' distribution, conductivity values at their locations and their amounts. Preserving only the Saint Venant's monopoles in respecting the "Venant's Condition" allows a significant improvement of the numerical results for the EEG forward problem.

## REFERENCES

1. Grech, R., T. Cassar, J. Muscat, et al., "Review on solving the inverse problem in EEG source analysis," *Jnl. of Neuroengineering and Rehabilitation*, Vol. 5, No. 1, 25, 2008.

2. Hämäläinen, M., R. Hari, R. J. Ilmoniemi, et al., "Magnetoencephalography — Theory, instrumentation, and applications to noninvasive studies of the working human brain," *Reviews of Modern Physics*, Vol. 65, No. 2, 413, 1993.
3. Tadel, F., S. Baillet, J. C. Mosher, et al., “Brainstorm: A user-friendly application for MEG/EEG analysis,” *Comput. Intelligence and Neuroscience*, Vol. 2011, 8, 2011.

4. Hallez, H., B. Vanrumste, R. Grech, et al., “Review on solving the forward problem in EEG analysis,” *Jnl. of Neuroengineering and Rehabilitation*, Vol. 4, No. 1, 46, 2007.

5. Wolters, C. H., H. Köstler, C. Möller, et al., “Numerical mathematics of the subtraction method for the modeling of a current dipole in EEG source reconstruction using finite element head models,” *SIAM Journal on Scientific Computing*, Vol. 30, No. 1, 24–45, 2007.

6. Bertrand, O., M. Thevenet, and F. Perrin, “3D finite element method in brain electrical activity studies,” *Biomagnetic Localization and 3D Modelling*, 154–171, 1991.

7. Shahid, S. and P. Wen, “Analytic and numeric evaluation of EEG forward problem using spherical volume conductor models,” *2010 IEEE/ICME International Conference on IEEE Complex Medical Engineering (CME)*, 28–33, 2010.

8. Wolters, C. H., A. Anwander, X. Tricoche, et al., “Influence of tissue conductivity anisotropy on EEG/MEG field and return current computation in a realistic head model: A simulation and visualization study using high-resolution finite element modeling,” *NeuroImage*, Vol. 30, No. 3, 813–826, 2006.

9. Zhang, Y., Z. Ren, and D. Lautru, “Finite element modeling of current dipoles using direct and subtraction methods for EEG forward problem,” *COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering*, Vol. 33, No. 1/2, 210, 2014.

10. Medani, T., D. Lautru, and Z. Ren, “Study of modeling of current dipoles in the finite element method for eeg forward problem,” *La 7ème Conférence Européenne sur les Méthodes Numériques en Electromagnétisme (NUMELEC 2012)*, Juillet, Marseille, France, 2012.

11. Zhang, Z., “A fast method to compute surface potentials generated by dipoles within multilayer anisotropic spheres,” *Physics in Medicine and Biology*, Vol. 40, No. 3, 335, 1995.

12. Buchner, H., Gr. Knoll, M. Fuchs, et al., “Inverse localization of electric dipole current sources in finite element models of the human head,” *Electroencephalography and Clinical Neurophysiology*, Vol. 102, No. 4, 267–278, 1997.

13. Yan, Y., P. Nunez, and R. T. Hart, “Finite-element model of the human head: Scalp potentials due to dipole sources,” *Medical and Biological Engineering and Computing*, Vol. 29, No. 5, 475–481, 1991.

14. Lew, S., C. H. Wolters, T. Dierkes, et al., “Accuracy and run-time comparison for different potential approaches and iterative solvers in finite element method based EEG source analysis,” *Applied Numerical Mathematics*, Vol. 59, No. 8, 1970–1988, 2009.

15. Drechsler, F., C. H. Wolters, T. Dierkes, et al., “A full subtraction approach for finite element method based source analysis using constrained Delaunay tetrahedralisation,” *NeuroImage*, Vol. 46, No. 4, 1055–1065, 2009.

16. Vorwerk, J., M. Clerc, M. Burger, and C. H. Wolters, “Comparison of boundary element and finite element approaches to the EEG forward problem,” *Biomedical Engineering*, Vol. 57(Suppl. 1), 795–798, 2012.

17. Vorwerk, J., J.-H. Cho, S. Rampp, H. Hamer, T. R. Knosche, and C. H. Wolters, “A guideline for head volume conductor modeling in EEG and MEG,” *NeuroImage*, Vol. 100, 590–607, 2014.

18. Fang, Q. and D. Boas, “Tetrahedral mesh generation from volumetric binary and gray-scale images,” *Proceedings of IEEE International Symposium on Biomedical Imaging 2009*, 1142–1145, 2009.

19. Liu, A. and B. Joe, “Relationship between tetrahedron shape measures,” *BIT Numerical Mathematics*, Vol. 34, No. 2, 268–287, 1994.

20. Medani, T., D. Lautru, Z. Ren, D. Schwartz, and G. Sou, “Modelling of brain sources using the modified Saint Venant’s method in FEM resolution of EEG forward problem,” *7th Annual International IEEE EMBS Conference on Neural Engineering*, France, Apr. 2015.

21. Gramfort, A., T. Papadopoulo, E. Olivi, and M. Clerc, “OpenMEEG: Opensource software for quasistatic bioelectromagnetics,” *BioMedical Engin.*, Vol. 45, No. 9, 2010.