The Discontinuous Galerkin Finite Element Method for Solving the MEG and the Combined MEG/EEG Forward Problem

In Electro- (EEG) and Magnetoencephalography (MEG), one important requirement of source reconstruction is the forward model. The continuous Galerkin finite element method (CG-FEM) has become one of the dominant approaches for solving the forward problem over the last decades. Recently, a discontinuous Galerkin FEM (DG-FEM) EEG forward approach has been proposed as an alternative to CG-FEM (Engwer et al., 2017). It was shown that DG-FEM preserves the property of conservation of charge and that it can, in certain situations such as the so-called skull leakages, be superior to the standard CG-FEM approach. In this paper, we developed, implemented, and evaluated two DG-FEM approaches for the MEG forward problem, namely a conservative and a non-conservative one. The subtraction approach was used as source model. The validation and evaluation work was done in statistical investigations in multi-layer homogeneous sphere models, where an analytic solution exists, and in a six-compartment realistically shaped head volume conductor model. In agreement with the theory, the conservative DG-FEM approach was found to be superior to the non-conservative DG-FEM implementation. This approach also showed convergence with increasing resolution of the hexahedral meshes. While in the EEG case, in presence of skull leakages, DG-FEM outperformed CG-FEM, in MEG, DG-FEM achieved similar numerical errors as the CG-FEM approach, i.e., skull leakages do not play a role for the MEG modality. In particular, for the finest mesh resolution of 1 mm sources with a distance of 1.59 mm from the brain-CSF surface, DG-FEM yielded mean topographical errors (relative difference measure, RDM%) of 1.5% and mean magnitude errors (MAG%) of 0.1% for the magnetic field. However, if the goal is a combined source analysis of EEG and MEG data, then it is highly desirable to employ the same forward model for both EEG and MEG data. Based on these results, we conclude that the newly presented conservative DG-FEM can at least complement and in some scenarios even outperform the established CG-FEM approaches in EEG or combined MEG/EEG source analysis scenarios, which motivates a further evaluation of DG-FEM for applications in bioelectromagnetism.

1. Introduction

Together with electroencephalography (EEG), magnetoencephalography (MEG) is a technique used to investigate brain activity. EEG and MEG are devoted to detect the electric potential distribution and the magnetic field generated by the brain, respectively, with a unique time resolution (Brette and Destexhe, 2012). An important topic in many applications of EEG and MEG is the source reconstruction, i.e., the identification of the sources in the brain responsible for the signals recorded at the head surface (EEG) or in a small distance from the head surface (MEG). Moreover, it has been shown (Fuchs et al., 1998; Aydin et al., 2015, 2017) that combined MEG/EEG employs the complementary information of both modalities providing source reconstructions that outperform the ones provided by each single modality. In order to compute MEG/EEG source reconstructions, i.e., to solve a related ill-posed inverse problem of MEG/EEG, the forward problem has to be solved. Since the accuracy of MEG/EEG inverse problem solutions depends strongly on the forward solution, it is fundamental to increase the accuracy of the latter (Brette and Destexhe, 2012). Furthermore, in a multi-modal MEG/EEG reconstruction it is desirable to use the same forward model for both EEG and MEG data. In the EEG case, the forward problem consists of the evaluation of the electric potential generated by a source located in the brain by solving an elliptic partial differential equation of second order (Wolters et al., 2007). In the MEG case, the magnetic field needs to be computed exploiting Biot-Savart's law, which depends on the EEG solution (Brette and Destexhe, 2012). In simplified scenarios, such as multi-layer sphere models with piecewise homogeneous conductivity, analytical solutions are available (Brette and Destexhe, 2012). In more realistic scenarios, e.g., realistically shaped head models, numerical methods have to be adopted. There is a large variety of numerical methods that can be employed, among them are boundary element methods (Mosher et al., 1999; Acar and Makeig, 2010; Gramfort et al., 2011; Stenroos and Sarvas, 2012), finite volume methods (Cook and Koles, 2006), finite difference methods (Wendel et al., 2008; Vatta et al., 2009; Montes-Restrepo et al., 2014) and finite element methods (FEMs) (Bertrand et al., 1991; Marin et al., 1998; Schimpf et al., 2002; Drechsler et al., 2009; Nüßing et al., 2016; Pursiainen et al., 2016). In this work, we deal with FEMs, which have shown high numerical accuracies with the possibility to model complex geometries and bioelectromagnetic properties (e.g., anisotropic conductivity) of the head. In more realistic simulations, an aspect that should be more carefully studied is the fulfilling of the conservation of charge law and its implications on the application at hand. For the EEG, this has been studied in Engwer et al. (2017), where it was shown that the phenomenon of skull leakages, which occur as a consequence of violating the conservation of charge law, can be overcome by using a discontinuous Galerkin FEM (DG-FEM) instead of a classical, continuous Galerkin FEM (CG-FEM). Leakage effects occur when a low conductive compartment of the head, i.e., the skull, is modeled too coarsely. This leads to scalp and cerebrospinal fluid elements being erroneously connected via single skull vertices or edges, a frequent case when segmenting, for example, children heads with thin skull compartments. Note that such skull leakage effects can also compromise the accuracy of transcranial electrical stimulation simulations (Miranda et al., 2006; Datta et al., 2013; Windhoff et al., 2013; Wagner et al., 2014), in a reciprocal sense (Wagner et al., 2016). As a further motivation, DG-FEM provides the basis for the implementation of a further improved method, the so-called unfitted discontinuous Galerkin FEM (Nüßing et al., 2016). This method combines the advantages of using a hexahedral mesh, whose generation pipeline is less complex than for the tetrahedral case, with a smooth representation of head tissue compartments.

In this work, we introduce the first application of DG-FEM for solving the MEG forward problem. As MEG solutions depend on EEG solutions, we implement a setup where the same method (CG- or DG-FEM) is adopted for both modalities, allowing for a combined EEG and MEG source reconstruction study. We analyze the accuracy of DG-FEM results in comparison to the CG-FEM results and investigate the propagation of the effects of conservation properties on the MEG results, as these effects have been proven to play an important role for EEG forward simulations in leaky scenarios. We will show that, in contrast to the EEG, the accuracy of the forward solution for the MEG is basically not affected by skull leakage effects. In fact, the accuracy of DG-FEM forward modeling is in the same range as for standard CG-FEM. However, because of the advantages on the EEG side, DG-FEM is an interesting new approach for combined MEG/EEG source reconstruction scenarios.

2. Theory

In this section, after a summary of the EEG and MEG background, the theory of Discontinuous Galerkin- (DG-) FEM for solving the MEG forward problem will be presented. As the MEG forward problem is build on the solution of the EEG forward problem, a brief section will be about the latter.

2.1. Background

Following Hämäläinen et al. (1993) and Brette and Destexhe (2012), the electric potential distribution and the resulting magnetic induction generated in the brain can be modeled through the quasistatic approximation of Maxwell's equations, when assuming that the permeability of the tissue in the head is that of the free space, i.e., μ = μ₀,

$\begin{array}{lll}\nabla \times \text{E}& =0,& (1\text{a})\end{array}$

$\begin{array}{lll}\nabla ·\text{E}& =\frac{\rho }{{ϵ}_0},& (1\text{b})\end{array}$

related to the electrical part, and

$\begin{array}{lll}\nabla \times \text{B}& ={\mu }_0\text{j},& (2\text{a})\end{array}$

$\begin{array}{lll}\nabla ·\text{B}& =0,& (2\text{b})\end{array}$

related to the magnetic part.

In Equation (2a) j represents the total current density produced by neuronal activity, which, in bioelectromagnetism (Hämäläinen et al., 1993; Brette and Destexhe, 2012), is split into two contributions,

$\begin{array}{ll}\text{j}(\text{r})={\text{j}}^p(\text{r})+{\text{j}}^s(\text{r}),& (3)\end{array}$

where j^p is the so called primary current, j^s the secondary or volume current and r ∈ ℝ³. In neuromagnetism, the primary current is widely represented as a mathematical point dipole (De Munck et al., 1988; Murakami and Okada, 2006),

$\begin{array}{ll}{\text{j}}^p(\text{r})=\text{M}·\delta (\text{r}-{\text{r}}_0),& (4)\end{array}$

where M ∈ ℝ³ stands for the dipolar moment and δ is the Dirac delta distribution, centered in the dipole position ${\mathbf{\text{r}}}_0\in {ℝ}^3$. The volume current is a passive current that is the result of the macroscopic electric field on charge carriers in the conducting medium (Hämäläinen et al., 1993; Brette and Destexhe, 2012), and

$\begin{array}{ll}{\text{j}}^s=\sigma \text{E}& (5)\end{array}$

holds true (Ohm's law), where σ indicates the conductivity profile of the conductive medium. While, for the mathematical point dipole, the primary current is present only at the source position, the secondary current flows passively everywhere in the medium.

2.1.1. The Forward Problem of EEG

To derive the EEG forward problem, we have to consider Equations (1a) and (2a). From Equation (1a) we deduce that there exists a potential u such that

$\begin{array}{ll}\text{E}=-\nabla u,& (6)\end{array}$

so that Equation (5) can be written as

$\begin{array}{ll}{\text{j}}^s=-\sigma \nabla u.& (7)\end{array}$

Applying the divergence to Equation (2a), we obtain

$\begin{array}{ll}\nabla ·\text{j}=0.& (8)\end{array}$

Combining Equations (3), (7), and (8), we get an inhomogeneous Poisson equation that, together with the homogeneous Neumann boundary condition, models the EEG forward problem:

$\begin{array}{lll}\nabla ·(\sigma \nabla u)=\nabla ·{\text{j}}^p,& \mathrm{\text{in}}\mathrm{ \Omega }\subseteq {ℝ}^3& (9)\end{array}$

$\begin{array}{lll}\sigma \nabla u·\text{n}=0,& \mathrm{\text{on}}\mathrm{ \partial }\Omega & (10)\end{array}$

where Ω is the volume conductor and n is the unit outer normal vector on ∂Ω.

2.1.2. Conservation Properties

A fundamental physical property of the EEG forward problem is the conservation of charge:

$\begin{array}{ll}{\displaystyle {\int }_{\partial K}{\text{j}}^s·\text{n}ds}={\displaystyle {\int }_{K}fdK,\forall K\subset \Omega },& (11)\end{array}$

where f = −∇ · j^p and K is a control volume in Ω.

For FEMs this property carries over to the discrete solution only if the test space contains the characteristic function, which is one in K and zero everywhere else. In general, a conforming discretization, like CG-FEM, does not guarantee this property, while the DG-FEM fulfills a discrete analog, see Remark 3.

2.1.3. The Forward Problem of MEG

The solution of the MEG forward problem consists in the computation of the magnetic induction (flux), Φ, generated by a dipolar source in the brain. The magnetic flux is computed from the magnetic field B (B-field):

$\begin{array}{ll}\Phi ={\displaystyle {\int }_S}\text{B}·d\text{s},& (12)\end{array}$

where S is the surface of the sensor.

Furthermore, following Biot-Savart's law, the B-field at a point r ∈ ℝ³ outside the domain Ω can be computed as

$\begin{array}{ll}\text{B}(\text{r})=\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_{\Omega }}\text{j}({\text{r}}^{\prime })\times \frac{\text{r}-{\text{r}}^{\prime }}{|\text{r}-{\text{r}}^{\prime }{|}^3}{d}^3{\text{r}}^{\prime },& (13)\end{array}$

(Hämäläinen et al., 1993; Brette and Destexhe, 2012).

When combining Equations (3), (13), and (4), one obtains (Hämäläinen et al., 1993; Brette and Destexhe, 2012):

$\begin{array}{l}\text{B}(\text{r})\stackrel{(3,13)}{=}\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_{\Omega }\left({\text{j}}^p(r^{\prime })+{\text{j}}^s(r^{\prime })\right)}\times \frac{\text{r}-r^{\prime }}{|\text{r}-r^{\prime }{|}^3}{d}^3{r}^{\prime }\\ \stackrel{(4)}{=}\frac{{\mu }_0}{4\pi }\text{M}\times \frac{\text{r}-{\text{r}}_0}{|\text{r}-{\text{r}}_0{|}^3}-\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_{\Omega }\sigma }\nabla u(r^{\prime })\times \frac{\text{r}-r^{\prime }}{|\text{r}-r^{\prime }{|}^3}{d}^3{r}^{\prime }\\ ={\text{B}}^p(\text{r})+{\text{B}}^s(\text{r}).& (14)\end{array}$

Namely, the B-field can be split into two contributions as well, the primary B-field B^p, which is calculated analytically for a mathematical point dipole in Equation (13), and the secondary B-field B^s, which has to be computed numerically when the electrical potential is computed numerically (since it depends on the electrical potential u inside the domain Ω).

2.1.4. The Forward Problem of MEG for Multi-layer Homogeneous Sphere Model

In simplified geometries, similarly to the EEG forward problem (Brette and Destexhe, 2012), there exist analytical solutions for the MEG forward problem (Sarvas, 1987; Ilmoniemi, 1995). Sarvas (1987) showed that the magnetic field outside a spherically symmetric conductor due to internal current sources does not depend on the profile of conductivity along the radius. He derived the following analytical MEG solution for a multi-layer homogeneous sphere model:

$\begin{array}{ll}\text{B}(\text{r})=\frac{{\mu }_0}{4\pi F^2}(F\text{M}\times {\text{r}}_0-\text{M}\times {\text{r}}_0·\text{r}\nabla F),& (15)\end{array}$

where a = r − r₀, a = |a|, r = |r|, $F=a(ra+r^2-\mathbf{\text{r}}·{\mathbf{\text{r}}}_0)$ and $\nabla F=(r^{-1}{a}^2+a^{-1}\mathbf{\text{a}}·\mathbf{\text{r}}+2a+2r)\mathbf{\text{r}}-(a+2r+a^{-1}\mathbf{\text{a}}·\mathbf{\text{r}}){\mathbf{\text{r}}}_0$.

Ilmoniemi (1995) could even demonstrate that radial anisotropy added to a spherically symmetric conductor does not affect the external magnetic field due to internal sources.

From Equation (15), one can deduce three important features of the analytical MEG solution for a multi-layer homogeneous sphere model and a point outside the model:

REMARK 1. Three main properties of analytical MEG solution for a multi-layer homogeneous sphere model and a measurement point outside the model (Sarvas, 1987):

1. the solution does not depend on the conductivity profile of the spherical model

2. if the source is radial, then the B-field outside Ω vanishes

3. the normal projection of the secondary component of the B-field gives a null contribution to the total B-field, i.e., B^s(r) · n = 0, for r outside Ω (Sarvas, 1987).

2.2. The EEG Forward Problem

In this section we will recall the concepts of CG- and DG-FEM for the EEG forward problem that are then needed in section 2.3 for the derivation of the two FEM based MEG forward approaches.

2.2.1. The Subtraction Approach

The mathematical point dipole model introduces a singularity on the right hand side of the PDE in Equation (9) that can be treated with the so-called subtraction approach (Bertrand et al., 1991; Awada et al., 1997; Marin et al., 1998; Wolters et al., 2007; Drechsler et al., 2009). The subtraction approach assumes that a non-empty neighborhood Ω^∞ around the source in r₀ can be found with homogeneous conductivity σ^∞. The conductivity tensor σ is then split into two parts,

$\begin{array}{ll}\sigma ={\sigma }^{\infty }+{\sigma }^{corr},& (16)\end{array}$

where σ^corr vanishes in Ω^∞. The potential u can also be split into two contributions,

$\begin{array}{ll}u=u^{\infty }+u^{corr}.& (17)\end{array}$

The so-called singularity potential u^∞ is the solution of the Poisson equation in an unbounded and homogeneous conductor with constant conductivity σ^∞, and it can be computed analytically. The correction potential u^corr becomes the unknown of a new Poisson equation:

$\begin{array}{lll}-\nabla ·(\sigma \nabla u^{corr})=\nabla ·({\sigma }^{corr}\nabla u^{\infty }),& \mathrm{\text{in}}\mathrm{ \Omega }\subseteq {ℝ}^3& (18)\end{array}$

$\begin{array}{lll}\sigma \nabla u^{corr}·\text{n}=-\sigma \nabla u^{\infty }·\text{n},& \mathrm{\text{on}}\mathrm{ \partial }\Omega & (19)\end{array}$

after embedding Equations (16) and (17) in Equations (9) and (10).

The conforming weak formulation of (18) and (19) presented in Wolters et al. (2007) reads: Find $u_h^{corr}\in W_h\subset H^1$ such that

$\begin{array}{l}{\displaystyle {\int }_{\Omega }}\sigma \nabla u_h^{corr}·\nabla v_{h}dx=-{\displaystyle {\int }_{\Omega }}{\sigma }^{corr}\nabla u^{\infty }·\nabla v_{h}dx\\ -{\displaystyle {\int }_{\partial \Omega }}{\sigma }^{\infty }\nabla u^{\infty }·\text{n}{v}_{h}ds& (20)\end{array}$

holds true, ∀v_h ∈ W_h. Choosing W_h as the space of piecewise linear, continuous functions give the classical CG-FEM.

The subtraction approach is theoretically well understood. A deep numerical analysis of the subtraction approach including proofs for uniqueness and existence has been carried out in Wolters et al. (2007) and Drechsler et al. (2009).

2.2.2. DG-FEM for the EEG Forward Problem

A discontinuous Galerkin- (DG-) FEM forward modeling approach has recently been proposed for the EEG by Engwer et al. (2017). In order to prepare our DG-FEM derivation for the MEG forward problem in section 2.3, we now recall some main properties of DG-FEM for the EEG. First, we introduce a volume triangulation ${{T}}_h(\Omega )$, which is a a finite collection of disjoint and open subsets forming a partition of Ω, where h corresponds to the mesh-width. Furthermore, the triangulation induces the internal skeleton

$\begin{array}{ll}{\Gamma }_{int}:=\left\{{\gamma }_{e,f}=\partial E_e\cap \partial E_f|E_e,E_f\in {{T}}_h(\Omega ),E_e\ne E_f,|{\gamma }_{e,f}|>0\right\}& (21)\end{array}$

and the skeleton Γ: = Γ_int ∪ ∂Ω. Let $V_h^l$ be the so-called broken polynomial space, that is defined as piecewise polynomial space on the partition ${{T}}_h(\Omega )$:

$\begin{array}{ll}{V}_h^l:=\left\{v\in L^2(\Omega ):v{|}_E\in P^l(E),\forall E\in {{T}}_h(\Omega )\right\},& (22)\end{array}$

where P^l denotes the space of polynomial functions of degree l ∈ ℕ. They describe functions that exhibit element-wise polynomial behavior but may be discontinuous across element interfaces. In the following we will assume that σ is constant on each element E_i and denote its value by σ_i.

Furthermore we recall the definition of jump of a function u on the intersection between two elements E_e and E_f of the triangulation ${{T}}_h(\Omega )$ with outer normal ${\mathbf{\text{n}}}_e\in {ℝ}^3$ and ${\mathbf{\text{n}}}_f\in {ℝ}^3$, respectively:

$\begin{array}{ll}⟦u⟧:=u{|}_{E_e}{\text{n}}_e+u{|}_{E_f}{\text{n}}_f\in {ℝ}^3.& (23)\end{array}$

Note that the normals n_e and n_f are opposing vectors, i.e., n_e = −n_f. In addition, the weighted average of u on the interface is defined as

$\begin{array}{ll}\left\{u\right\}:=\frac{{\sigma }_f}{{\sigma }_e+{\sigma }_f}u{|}_{E_e}+\frac{{\sigma }_e}{{\sigma }_e+{\sigma }_f}u{|}_{E_f}.& (24)\end{array}$

The DG-FEM for solving Equations (18) and (19) then reads: Find $u_h^{corr}\in V_h^l$ such that

$\begin{array}{ll}a(u_h^{corr},v_h)+J(u_h^{corr},v_h)=l(v_h),\forall v_h\in V_h^l,& (25)\end{array}$

with

$\begin{array}{l}a(u_h^{corr},v_h)={\displaystyle {\int }_{\Omega }}\sigma \nabla u_h^{corr}·\nabla v_{h}dx-{\displaystyle {\int }_{{\Gamma }_{int}}}\left\{\sigma \nabla u_h^{corr}\right\}·⟦v_h⟧ds\\ -{\displaystyle {\int }_{{\Gamma }_{int}}}\left\{\sigma \nabla v_h\right\}·⟦u_h^{corr}⟧ds,\\ J(u_h^{corr},v_h)=\eta {\displaystyle {\int }_{{\Gamma }_{int}}}\frac{{\widehat{\sigma }}_{\gamma }}{h_{\gamma }}⟦u_h^{corr}⟧·⟦v_h⟧ds,\end{array}$

and

$\begin{array}{l}l(v_h)=-{\displaystyle {\int }_{\Omega }}{\sigma }^{corr}\nabla u^{\infty }·\nabla v_{h}dx-{\displaystyle {\int }_{\partial \Omega }}{\sigma }^{\infty }\nabla u^{\infty }·\text{n}{v}_{h}ds\\ +{\displaystyle {\int }_{{\Gamma }_{int}}}\left\{{\sigma }^{corr}\nabla u^{\infty }\right\}·⟦v_h⟧ds,& (26)\end{array}$

where h_γ and ${\widehat{\sigma }}_{\gamma }$ denote local definitions of the mesh width and the electric conductivity on an edge γ, respectively; while η is a penalty parameter.

If Equation (25) has been solved toward the correction potential $u_h^{corr}$, then the full potential u_h can be computed as $u_h=u_h^{corr}+u^{\infty }$.

REMARK 2. (Discrete Properties) The proposed discretization Equation (25) is consistent and adjoint-consistent with the strong problem Equations (18) and (19), and for a sufficiently large constant η > 0 it has a unique solution.

REMARK 3. (Conservation Property) For any control volume $K\in {{T}}_h(\Omega )$, Equation (25) fulfills a discrete conservation property

$\begin{array}{ll}{\displaystyle {\int }_{\partial K}}{\text{j}}_h^{corr,DG}·\text{n}ds={\displaystyle {\int }_K}{f}^{corr}dx,& (27)\end{array}$

with the discrete electric flux ${\mathit{\text{j}}}_h^{corr,DG}=\left\{\sigma \nabla u_h^{corr}\right\}-\eta \frac{{\widehat{\sigma }}_{\gamma }}{h_{\gamma }}⟦u_h^{corr}⟧$ and f^corr = −∇ · σ^corr∇u^∞. For h → 0, the jump $⟦u_h^{corr}⟧$ vanishes and the discrete flux ${\mathit{\text{j}}}_h^{corr,DG}$ converges to the flux j^corr = σ∇u^corr.

As we will see in more detail in section 2.3.3, for the MEG problem the main quantity of interest is the electric flux, needed to compute the B-field. This flux, again, is closely related to the conservation of charge property.

2.3. The MEG Forward Problem

In this section, CG- and DG-FEM formulations for the MEG forward problem are derived, when a conservative or a non-conservative flux expression is adopted along with the subtraction approach.

2.3.1. The Subtraction Approach in the MEG Case

In this section, we will focus on the expression of the secondary B-field B^s, as the primary B-field B^p is analytically computable.

When inserting Equation (17) into Biot-Savart's law Equation (14), we obtain the following expression for the secondary B-field:

$\begin{array}{l}{\text{B}}^s(\text{r})=-\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_{\Omega }}(\sigma \nabla (u^{\infty }+u^{corr})({\text{r}}^{\prime }))\times \frac{\text{r}-{\text{r}}^{\prime }}{|\text{r}-{\text{r}}^{\prime }{|}^3}{d}^3{\text{r}}^{\prime }\\ =-\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_{\Omega }}\sigma \nabla u^{\infty }({\text{r}}^{\prime })\times \frac{\text{r}-{\text{r}}^{\prime }}{|\text{r}-{\text{r}}^{\prime }{|}^3}{d}^3{\text{r}}^{\prime }\\ -\frac{{\mu }_0}{4\pi }{\displaystyle {\int }_{\Omega }}\sigma \nabla u^{corr}({\text{r}}^{\prime })\times \frac{\text{r}-{\text{r}}^{\prime }}{|\text{r}-{\text{r}}^{\prime }{|}^3}{d}^3{\text{r}}^{\prime }\\ ={\text{B}}_{\infty }^s(\text{r})+{\text{B}}_{corr}^s(\text{r}).& (28)\end{array}$

Both ${\mathbf{\text{B}}}_{\infty }^s$ and ${\mathbf{\text{B}}}_{corr}^s$ are computed by numerical integration on the volume Ω.

Since we are dealing with numerical integration, both u^∞ and u^corr are projected in a discrete space, W_h, i.e.,

$\begin{array}{ll}{u}_h^{\infty }={\displaystyle \sum _i}{u}_i^{\infty }{\phi }_i,& (29)\end{array}$

and

$\begin{array}{ll}{u}_h^{corr}={\displaystyle \sum _i}{u}_i^{corr}{\phi }_i,& (30)\end{array}$

where (φ_i)_i represent a basis of the discrete space W_h, while $u_h^{\infty }$ and $u_h^{corr}$ are the discrete representations of u^∞ and u^corr, respectively. Note that u^∞ has an analytical expression. The description of the discretization process in both the CG- and DG-FEM schemes is the content of the following sections.

2.3.2. CG-FEM MEG Forward Problem

In a CG-FEM approach the following expression of the electric flux (${\mathbf{\text{j}}}_h^{corr,CG}$) is considered:

$\begin{array}{l}{\text{j}}_h^{corr,CG}=\sigma \nabla u_h^{corr}\\ ={\displaystyle \sum _i}{u}_i^{corr}\nabla {\phi }_i.& (31)\end{array}$

where (φ_i)_i is a collection of hat functions, basis of W_h.

The discretization of ${\mathbf{\text{B}}}_{\infty }^s$ and ${\mathbf{\text{B}}}_{corr}^s$ (i.e., ${\mathbf{\text{B}}}_{\infty ,h}^s$ and ${\mathbf{\text{B}}}_{corr,h}^s$) are then:

$\begin{array}{ll}{\text{B}}_{\infty ,h}^s(\text{r})=-\frac{{\mu }_0}{4\pi }{\displaystyle \sum _i}{u}_i^{\infty }{\displaystyle {\int }_{\Omega }}\sigma \nabla {\phi }_i({\text{r}}^{\prime })\times \frac{\text{r}-{\text{r}}^{\prime }}{|\text{r}-{\text{r}}^{\prime }{|}^3}{d}^3{\text{r}}^{\prime },& (32)\end{array}$

and

$\begin{array}{ll}{\text{B}}_{corr,h}^s(\text{r})=-\frac{{\mu }_0}{4\pi }{\displaystyle \sum _i}{u}_i^{corr}{\displaystyle {\int }_{\Omega }}\sigma \nabla {\phi }_i({\text{r}}^{\prime })\times \frac{\text{r}-{\text{r}}^{\prime }}{|\text{r}-{\text{r}}^{\prime }{|}^3}{d}^3{\text{r}}^{\prime },& (33)\end{array}$

respectively. Note that ${(u_i^{corr})}_i$ are given from the EEG forward computation.

If we call c_n the center of the nth coil, then the discrete ${\mathbf{\text{B}}}_{\infty ,h}^s$ and ${\mathbf{\text{B}}}_{corr,h}^s$ evaluated in c_n are:

$\begin{array}{ll}{\text{B}}_{\infty ,h}^s({\text{c}}_n)={\displaystyle \sum _i}{u}_i^{\infty }{\displaystyle \underset{:\; =S_{ni}^{\infty }}{{\displaystyle \underbrace{\left(-\frac{{\mu }_0}{4\pi }\right){\displaystyle {\int }_{\Omega }}\sigma \nabla {\phi }_i({\text{r}}^{\prime })\times \frac{{\text{c}}_n-{\text{r}}^{\prime }}{|{\text{c}}_n-{\text{r}}^{\prime }{|}^3}{d}^3{\text{r}}^{\prime }}}}},& (34)\end{array}$

and

$\begin{array}{ll}{\text{B}}_{corr,h}^s({\text{c}}_n)={\displaystyle \sum _i}{u}_i^{corr}{\displaystyle \underset{:\; =S_{ni}^{corr,CG}}{{\displaystyle \underbrace{\left(-\frac{{\mu }_0}{4\pi }\right){\displaystyle {\int }_{\Omega }}\sigma \nabla {\phi }_i({\text{r}}^{\prime })\times \frac{{\text{c}}_n-{\text{r}}^{\prime }}{|{\text{c}}_n-{\text{r}}^{\prime }{|}^3}{d}^3{\text{r}}^{\prime }}}}},& (35)\end{array}$

respectively.

$S^{\infty }={(S_{ni}^{\infty })}_{n,i}$ and $S^{corr,CG}={(S_{ni}^{corr,CG})}_{n,i}$ are the secondary magnetic field integration matrices. Equations (34) and (35) can be rewritten into matrix equations,

$\begin{array}{ll}{\text{B}}_{\infty ,h}^s=S^{\infty }{\text{u}}^{\infty },& (36)\end{array}$

and

$\begin{array}{ll}{\text{B}}_{corr,h}^s=S^{corr,CG}{\text{u}}^{corr},& (37)\end{array}$

respectively.

An alternative treatment of ${\mathbf{\text{B}}}_{corr}^s$ involves the already mentioned conservation of charge property, whose fulfillment is not guaranteed for the flux in Equation (31). Section 2.3.3 is dedicated to the description of this alternative.

2.3.3. Conservative Flux DG-FEM MEG Forward Problem

Following Equation (31), we can consider the analogous formula for the electric flux in the DG-FEM scheme, i.e.,

$\begin{array}{ll}\sigma \nabla u_h^{corr}={\displaystyle \sum _i}{u}_i^{corr}\nabla {\phi }_i,& (38)\end{array}$

where (φ_i)_i is a basis of $V_h^1$.

As already mentioned, in general this discrete formulation of the flux does not verify the conservation of charge property. Conversely and despite the CG-FEM case, in the DG-FEM approach we can consider another expression of the discrete electric flux, i.e.,

$\begin{array}{ll}{\text{j}}_h^{corr,DG}=\left\{\sigma \nabla u_h^{corr}\right\}-\eta \frac{{\widehat{\sigma }}_{\gamma }}{h_{\gamma }}⟦u_h^{corr}⟧& (39)\end{array}$

that verifies the conservation of charge law, as described in Remark 3.

The main idea is to embed this conservative current (or conservative flux, i.e., flux fulfilling the property of conservation of charge) in the computation of ${\mathbf{\text{B}}}_{corr}^s$. We have to notice that ${\mathbf{\text{j}}}_h^{corr,DG}$ is defined only on the internal skeleton Γ_int (Equation 21) and not in the entire volume Ω. In order to integrate ${\mathbf{\text{j}}}_h^{corr,DG}$ when computing ${\mathbf{\text{B}}}_{corr}^s$ (Equation 28), we need to project the current into the volume. One way to do so is to interpolate ${\mathbf{\text{j}}}_h^{corr,DG}$ in the space of the lowest-order Raviart Thomas function (RT₀). RT₀ is H(div)- conforming and its degrees of freedom (DOFs) are the evaluations of the basis functions along the projections normal to the faces of each element, exactly where ${\mathbf{\text{j}}}_h^{corr,DG}$ is defined. The space H(div; Ω) is defined as:

$\begin{array}{ll}H(div;\Omega ):=\left\{\text{v}\in L^2{(\Omega )}^3:\nabla ·\text{v}\in L^2(\Omega )\right\},& (40)\end{array}$

and RT₀ as (Nédélec, 1980; Fortin and Brezzi, 1991):

$\begin{array}{ll}R{T}_0({{T}}_h(\Omega )):\text{}=\{\text{v}\in H(div;\Omega ):(\nabla ·\text{v}){|}_E\in P^0(E),\forall E\in {{T}}_h(\Omega )\}.& (41)\end{array}$

As we are considering hexahedral elements, P^l(E) in Equation (22) is ℚ^l(E), defined as:

$\begin{array}{ll}{\mathbb{Q}}^l(E)=\mathrm{\text{span}}\left\{{\Pi }_{i=1}^3{x}_i^{{\alpha }_i}:\text{x}\in E,\alpha \in {\mathbb{N}}^3,\mathrm{\text{max }}{\alpha }_i\le l\right\},& (42)\end{array}$

therefore also in Equation (41), we have P⁰ = ℚ⁰. For a regular, hexahedral mesh with edge length h, as in our case, a RT₀ basis function ψ_k is supported on the two hexahedral elements $E_e,E_f\in {{T}}_h(\Omega )$ sharing the face f_k = E_e ∩ E_f with normal vector n_k and centroid ${\overline{\mathbf{\text{x}}}}_k$. It can be defined by

$\begin{array}{ll}{\mathit{\psi }}_k(\text{x})=\{\begin{array}{ll}\left(1+\frac{(\text{x}-{\overline{\text{x}}}_k)·{\text{n}}_k}{h}\right){\text{n}}_k,\hfill & \text{if }\text{x}\in {\overline{E}}_e\cap {\overline{E}}_f\hfill \\ 0,\hfill & \text{otherwise}\text{.}\hfill \end{array}& (43)\end{array}$

For more insights see Fortin and Brezzi (1991) and Nédélec (1980), and Figure 1, where the basis function ψ_k has been visualized.

FIGURE 1

Figure 1. Visualization of a zeroth-order Raviart-Thomas basis function (Right) and its support (Left). The support is made of two hexahedral elements E_e and E_f, which are sharing the face f_k with unit outer normal n_k. The vector valued function is equal to 1 · n_k on the face f_k and it decays when reaching the other parallel faces.

For the discretization of ${\mathbf{\text{B}}}_{corr}^s$ we can start from observing that the conservative flux $\Psi (u_h^{corr})$ is a function of $L^2({\Gamma }_{int})$ which depends on the potential $u_h^{corr}$:

$\begin{array}{ll}\text{Ψ}(u_h^{corr})={\text{j}}_h^{corr,DG}=\left\{\sigma \nabla u_h^{corr}\right\}-\eta \frac{{\widehat{\sigma }}_{\gamma }}{h_{\gamma }}⟦u_h^{corr}⟧\in L^2({\Gamma }_{int}).& (44)\end{array}$

If (φ_i)_i is a basis of $V_h^1$, then the correction potential can be written as

$\begin{array}{ll}{u}_h^{corr}={\displaystyle \sum _i}{u}_i^{corr}{\phi }_i,& (45)\end{array}$

and, due to linearity, we have

$\begin{array}{ll}\text{Ψ}(u_h^{corr})={\displaystyle \sum _i}{u}_i^{corr}\text{Ψ}({\phi }_i).& (46)\end{array}$

If we now apply the projection Π_RT0 into RT₀ to $\Psi (u_h^{corr})$ and we exploit again linearity, we obtain:

$\begin{array}{ll}{\Pi }_{RT0}(\text{Ψ}(u_h^{corr}))={\displaystyle \sum _i}{u}_i^{corr}{\Pi }_{RT0}(\text{Ψ}({\phi }_i))\in L^2(\Omega ).& (47)\end{array}$

Finally, ${\mathbf{\text{B}}}_{corr}^s$ can be then approximated as follows:

$\begin{array}{ll}{\text{B}}_{corr}^s(\text{r})\simeq -\frac{{\mu }_0}{4\pi }{\displaystyle \sum _i}{u}_i^{corr}{\displaystyle {\int }_{\Omega }}{\Pi }_{RT0}(\text{Ψ}({\phi }_i))({\text{r}}^{\prime })\times \frac{\text{r}-{\text{r}}^{\prime }}{|\text{r}-{\text{r}}^{\prime }{|}^3}{d}^3{\text{r}}^{\prime }.& (48)\end{array}$

If we call c_n the center of the nth coil, then the discretization of ${\mathbf{\text{B}}}_{corr}^s$ evaluated in c_n reads,

$\begin{array}{ll}{\text{B}}_{corr}^s({\text{c}}_n)\simeq {\displaystyle \sum _i}{u}_i^{corr}{\displaystyle \underset{:=S_{ni}^{corr,DG}}{{\displaystyle \underbrace{\left(-\frac{{\mu }_0}{4\pi }\right){\displaystyle \underset{\Omega }{\int }}{\Pi }_{RT0}(\text{Ψ}({\phi }_i))({\text{r}}^{\prime })\times \frac{{\text{c}}_n-{\text{r}}^{\prime }}{|{\text{c}}_n-{\text{r}}^{\prime }{|}^3}{d}^3{\text{r}}^{\prime }}}}}.& (49)\end{array}$

$S^{corr,DG}={(S_{ni}^{corr,DG})}_{n,i}$ is the secondary magnetic field integration matrix related to the DG-FEM scheme. Equation (49) can be rewritten into a matrix equation,

$\begin{array}{ll}{\text{B}}_{corr,h}^s=S^{corr,DG}{\text{u}}^{corr},& (50)\end{array}$

where ${\mathbf{\text{B}}}_{corr,h}^s$ represents the discretization of ${\mathbf{\text{B}}}_{corr}^s$.

REMARK 4. The projection of the i^th basis function of the space $V_h^1$ can be described as:

$\begin{array}{ll}{\Pi }_{R{T}_0}(\text{Ψ}({\phi }_i))={\displaystyle \sum _k}{\alpha }_i^k{\text{ψ}}_k,& (51)\end{array}$

where (ψ_k)_k form a basis of RT₀ and ${\alpha }_i^k$ are the DOFs, which can be derived as

$\begin{array}{ll}{\alpha }_i^k=\text{Ψ}({\phi }_i({\bar{\text{x}}}_k))·{\text{n}}_k,& (52)\end{array}$

with ${\overline{\mathbf{\text{x}}}}_k$ and n_k the centroid and the external normal of the face f_k, respectively (see Figure 1).

2.4. Transfer Matrix Approach

As described in the next section, MEG forward computations will be carried out for a large number of dipole sources. In order to speed up the many numerically expensive computations of the secondary B-field B^s for all of these sources, following Wolters et al. (2004), we adapted and implemented transfer matrix approaches for all three presented FEM-based MEG forward modeling schemes.

If Ku = j represents the resulting linear system of the EEG forward computation discretization, we can formally write

$\begin{array}{ll}\text{u}=K^{-1}\text{j}.& (53)\end{array}$

If we combine Equations (53) and (50), we obtain

$\begin{array}{ll}{\text{B}}_{corr,h}^s=S\text{u}=S(K^{-1}\text{j})=(S{K}^{-1})\text{j}=B_{MEG}\text{j},& (54)\end{array}$

where S is a generic secondary magnetic field integration matrix. B_MEG is the so-called MEG transfer matrix and allows computing ${\mathbf{\text{B}}}_{corr,h}^s$ with a matrix-vector multiplication, instead of solving the EEG forward problem and applying S.

To compute $B_{MEG}=S{K}^{-1}$, we can multiply its definition by K from the right and then transpose it.

Using the symmetry of K, we arrive at the following matrix equation,

$\begin{array}{ll}K{B}_{MEG}^t=S^t,& (55)\end{array}$

which can be solved for each row of S (column of S^t).

3. Methods

3.1. Implementation

We implemented the CG-FEM and the two DG-FEM approaches [non-conservative Equation 38 and conservative flux (Equation 39)] for the MEG forward problem in the Distributed and Unified Numerics Environment (DUNE)¹ (Blatt and Bastian, 2007; Bastian et al., 2008a,b). DUNE is a modular open source C++ library for solving partial differential equations with mesh-based methods. In particular, we used the DUNE-ALUGrid module (Alkämper et al., 2016) for the representation of hexahedral meshes and the DUNE-PDELab module (Bastian et al., 2010) for the discretization of the partial differential equations. All the newest implementations have been gathered in the module called DUNEuro², a special module dedicated to solve PDEs in neuroscience.

3.2. Volume Conductor Models

For numerical accuracy tests of our new CG- and DG-FEM implementations, we generated 4-layer homogeneous sphere models for which an analytical solution for the MEG exists (see section 2.1.4). We used four compartments with different conductivities in order to evaluate if, besides the analytical solution in Equation (15), also our numerical implementations show conductivity-independence of MEG in spherical volume conductors and because the four compartment model is closer to a realistic head model as shown in Figure 10. The four compartments, whose radii and conductivities are shown in Table 1 (same parametrization as in Engwer et al., 2017), are rough approximations for skin, skull, cerebrospinal fluid (CSF) and brain compartments. The spherical domain is represented via hexahedral meshes with three different resolutions, namely 4, 2, and 1 mm. In this work, we focused on hexahedral meshes in order to study the scenario where the combination of thin skull structures and insufficient hexahedral mesh resolutions might result in so-called skull leakages, as in Engwer et al. (2017). The numbers of vertices and elements of these meshes are shown in Table 2.

TABLE 1

Table 1. Four compartment sphere model.

TABLE 2

Table 2. Parameters (from left to right) of the regular hexahedral meshes of the 4-layer sphere models used for validation purposes: segmentation resolution (Segm. Res.), mesh width (h), number of vertices and number of elements.

3.3. Sources and Sensors

As only tangential orientation components produce an MEG signal in a multi-layer sphere model (section 2.1.4), we generated 8,000 dipoles with purely tangential orientations and unit strengths. The sources were uniformly distributed inside the brain compartment on spherical surfaces with 8 different logarithmically scaled eccentricities reported in Table 3. A source with eccentricity value of 0 is positioned in the center of the sphere, while a source with eccentricity value of 1 belongs to the surface separating brain and CSF compartments. The logarithmic scaling was chosen, since it is well known that numerical errors of the subtraction approach increase with decreasing distance of a source to the next conductivity jump (Wolters et al., 2007; Drechsler et al., 2009). We therefore expect larger numerical errors especially for the sources at the highest eccentricity of 0.9873, which only have a distance of 0.99 mm to the CSF compartment.

TABLE 3

Table 3. Source eccentricities and corresponding distances to the CSF compartment.

As the cortex has a thickness of 4 to 2 mm (Hämäläinen et al., 1993; Murakami and Okada, 2006) and the sources are located in the center of the gray matter, the sources which are most important to analyze are those with a distance of 2 to 1 mm to the CSF compartment. Therefore we will focus on the results of sources whose eccentricities are between 0.9642 (2.79 mm from the CSF compartment) and 0.9873 (0.99 mm from the CSF compartment) and especially on those with the middle value of this range, i.e., 0.9796 (1.59 mm from the CSF compartment). Furthermore, in praxis (and for the realistic head model used in this study, section 3.6), sources are usually placed so that at least one layer of elements is between the source element and the conductivity jump, which is fulfilled for the considered eccentricities ⩽ 0.9873 in the 1 mm model (seg_1_res_1) and ⩽ 0.9642 in the 2 mm model (seg_2_res_2). See Table 3 for details on the eccentricities and the corresponding distance from the CSF compartment.

With regard to the MEG sensors, we used 256 point-magnetometers outside the sphere model at a fixed radius of 110 mm (see Figure 2).

FIGURE 2

Figure 2. Visualization of the 256 point-magnetometers used in the sphere model analysis. Radially (Left) and tangentially (Middle) oriented point-magnetometers have been employed exclusively in section 4.1.1, while in all other studies all the three Cartesian components (Right) of the vector fields B^p, B^s, and B have been considered.

3.4. Error Measures

We will use the two error metrics that are commonly used for validating EEG and MEG forward approaches (Meijs et al., 1989; Bertrand et al., 1991; Marin et al., 1998; van den Broek et al., 1998; Schimpf et al., 2002; Wolters et al., 2007), namely, the relative difference measure (RDM%) for topographical errors:

$\begin{array}{ll}RDM\%({\text{f}}_{ana},{\text{f}}_{num})=50{\Vert \frac{{\text{f}}_{num}}{||{\text{f}}_{num}|{|}_2}-\frac{{\text{f}}_{ana}}{||{\text{f}}_{ana}|{|}_2}\Vert }_2,& (56)\end{array}$

and magnitude error (MAG%):

$\begin{array}{ll}MAG\%({\text{f}}_{ana},{\text{f}}_{num})=100\left(\frac{||{\text{f}}_{num}|{|}_2}{||{\text{f}}_{ana}|{|}_2}-1\right),& (57)\end{array}$

where f is either the secondary B-field B^s or the full B-field B. Note that we considered vector-magnetic fields (B^p, B^s, B) without projecting them into radial nor tangential directions, i.e., without distinguishing between radial and tangential point-magnetometers.

Statistical results of numerical accuracies will be visualized with mean curves and boxplots (see Figures 5–9). In the boxplots, the analysis includes maximum and minimum, indicated by upper and lower error bars, and thereby the total range (TR). Furthermore, it includes the interval between upper and lower quartile, i.e., the interquartile range (IQR), which is marked by a box with a black dash showing the median.

3.5. A Leaky Model

The study of Engwer et al. (2017) showed for EEG forward scenarios that DG-FEM can considerably outperform CG-FEM in skull leakage models, where the sphere model is discretized with a hexahedral mesh of 2 mm resolution and where at the same time the thickness of the skull compartment is deliberately reduced down to 2 mm (seg_2_res_2_r82), so that it contains many leaky points, i.e., vertices belonging to both an element labeled as skin and an element labeled as CSF or brain. See Table 4 for details. Note that real skull holes are not investigated in this study. The sources are the same as previously described. Even when not expecting similarly substantial error reductions on the MEG side, we used here the same leakage models as in Engwer et al. (2017) to investigate the influence of skull leakages on the presented CG- and DG-FEM MEG approaches.

TABLE 4

Table 4. Parameters (from left to right) of the regular hexahedral meshes of the 4-layer sphere models used to investigate the influence of skull leakages on the presented CG- and DG-FEM MEG approaches: segmentation resolution (Segm. Res.), mesh width (h), outer radius of the skull (mm) and number of leaky points.

3.6. A Realistic Head Model

As a proof of concept, we computed one MEG forward solution using the DG-FEM approach in a more realistic scenario. Based on MRI recordings of a human head, a segmentation considering six tissue compartments (white matter, gray matter, cerebrospinal fluid, skull compacta, skull spongiosa, and skin) that includes realistic skull openings such as the foramen magnum and the optic nerve canal was generated. Based on this segmentation, a six-compartment realistically shaped head model was built, a hexahedral mesh of 2 mm resolution resulting in 508,412 vertices and 484,532 elements (Figure 10). As this model was not corrected for leakages, 1,164 vertices belonging to both CSF and skin elements were found. These leaky points were mainly located at the temporal bone. More details about the model and its generation process can be found in Engwer et al. (2017). Locations and orientations of the sensors were chosen accordingly to the CTF machine (OMEGA2005, CTF, VSM MedTech Ltd., Canada), see Figure 10.

4. Results

In this section, the results relative to the evaluation and validation in multi-layer homogeneous sphere will be presented, followed by the results of one forward computation on a realistically shaped head model.

4.1. Validations and Evaluations in Spherical Volume Conductor Models

In this section, we will validate, compare and evaluate the three developed and implemented approaches for the MEG forward problem, namely the CG-FEM and the DG-FEM with non-conservative (Equations 31, 38) and conservative flux (Equation 39), in spherical volume conductor models.

4.1.1. Preparatory Work Using Analytical Approach

To recall the most important symmetry properties of the MEG forward problem in spherical volume conductor models, to prepare the numerical studies below and to enable an easier interpretation of their results, we first tested and visualized the properties of the MEG analytical solution for a multi-layer homogeneous sphere model, as reported in Remark 1. Here, we consider radial and tangential point-magnetometers, i.e., we have projected the B-field (B^p,B^s,B) onto the radial n and tangential t directions at sensor locations (Figure 2).

In Figure 3A, we compared, for the tangentially-oriented sources at logarithmically scaled eccentricities and the 256 radial point-magnetometers, the L² norm of primary B^p (in pink) and secondary B^s (in blue) B-fields, i.e.,

$\begin{array}{ll}||{\text{B}}^p·\text{n}|{|}_2,||{\text{B}}^s·\text{n}|{|}_2.& (58)\end{array}$

We notice that the only contribution to radial point-magnetometers is given by the primary component of the B-field, B^s, as proven in Sarvas (1987).

FIGURE 3

Figure 3. (A) Analytical solutions in spherical volume conductor model for radial point-magnetometers: L² norm of the primary (B^p, pink) and secondary (B^s, blue) B-fields (see Equation 58) for tangentially-oriented sources at logarithmically scaled eccentricities. Values are expressed in Tesla T. (B) Analytical solutions in spherical volume conductor model for radial point-magnetometers: L² norm of the radial full B-field component relative to the one for the most eccentric source (see Equation 59) for tangentially-oriented sources at logarithmically scaled eccentricities.

In Figure 3B, we plotted the L² norm of the full B-field for radial point-magnetometers normalized to the maximum over all tested sources, which is achieved for the most eccentric source, i.e.,

$\begin{array}{ll}\frac{||\text{B}·\text{n}|{|}_2}{\mathrm{\text{max}}||\text{B}·\text{n}|{|}_2}.& (59)\end{array}$

We can see how the magnitude of the full B-field increases for sources with an increasing eccentricity.

In Figure 4, we investigated the analytical solutions in the spherical volume conductor model for tangential point-magnetometers (Figure 2, middle). The Figure shows the L² norm of the primary (in pink) and secondary (in blue) tangential B-field components, i.e.,

$\begin{array}{ll}||{\text{B}}^p·\text{t}|{|}_2,||{\text{B}}^s·\text{t}|{|}_2,& (60)\end{array}$

for tangentially-oriented sources at different eccentricities.

FIGURE 4

Figure 4. Analytical solutions in spherical volume conductor model for tangential point-magnetometers: L² norm of the primary (B^p, pink) and secondary (B^s, blue) B-fields (see Equation 60) for tangentially-oriented sources at logarithmically scaled eccentricities. Values are expressed in Tesla T.

In this Figure we can see that, for tangential point-magnetometers, the deeper the sources are, the more the primary and secondary B-fields give identical contributions, but with opposite signs, to the full B-field, i.e., they more and more cancel each other out. Toward the sphere center, sources become more and more radial and the full B-field goes down to zero. However, as Figure 4 also shows, with increasing source eccentricity the relative contribution of the primary tangential B-field component increases when compared to the secondary B-field component. The tangential full B-field projection (i.e., B · t) and, together with it, the difference between primary and secondary tangential B-field components (i.e., B^p · t and B^s · t) thus increase with increasing source eccentricity.

4.1.2. FEM Study 1: Conservative vs. Non-conservative Flux for DG Approach

We now turn our interest to the validation and evaluation of our implemented new numerical FEM approaches for the MEG forward problem in spherical models. We will only consider tangentially-oriented sources for the validations and evaluations in the next sections, because, as seen in section 4.1.1, radial sources do not produce any magnetic field outside spherical volume conductor models. Following Equations (28) and (14), we will from now on measure errors of the vector fields B^s (Figures 5, 6, 7, 9) and B (Figure 8). These errors thus include parts from the radial and the two tangential sensor orientations and thus enable an overall view on the MEG forward modeling accuracy. On the one hand, radially-oriented sensor orientations are dominant in realistic MEG sensor configurations (see Figure 10), while on the other hand, and as seen in section 4.1.1, because of the cancellation effect of primary and secondary B-fields, tangentially-oriented sensor orientations are especially delicate numerical test-cases.

FIGURE 5

Figure 5. Accuracy comparison for secondary B-field B^s computation (Equation 28) between DG-FEM with non-conservative flux (Equation 38, in red) and DG-FEM with the conservative flux (Equation 39, in green) in a 4 mm hexahedral sphere model: visualized are the means (Left column) and the boxplots (Right column) of the RDM% (Top row) and MAG% (Bottom row), for tangentially oriented sources at logarithmically-scaled eccentricities. Dipoles not belonging to the brain compartment are excluded from the statistics. The dashed green line represents the eccentricity of 4 mm distance to the brain-CSF boundary. Note the different scaling of the y-axes (Top row).

FIGURE 6

Figure 6. Validation and convergence analysis for secondary B-field B^s computation (Equation 28) of DG-FEM with conservative flux (Equation 39) in a 4 mm (green), 2 mm (red) and 1 mm (blue) hexahedral sphere model: visualized are the means (Left column) and the boxplots (Right column) of the RDM% (Top row) and MAG% (Bottom row), for tangentially oriented sources at logarithmically-scaled eccentricities. Dipoles not belonging to the brain compartment are excluded from the statistics. Dashed lines represent the eccentricities of 4 mm (green), 2 mm (red) and 1 mm (blue) distances to the brain-CSF boundary. Note the different scaling of the y-axes (Top row).

FIGURE 7

Figure 7. Accuracy comparison for secondary B-field B^s computation (Equation 28) between CG-FEM (in warm colors) and DG-FEM with the conservative flux (in cold colors), for different mesh resolutions: visualized are the means (Left column) and the boxplots (Right column) of the RDM% (Top row) and MAG% (Bottom row), for tangentially oriented sources at logarithmically-scaled eccentricities. Dipoles not belonging to the brain compartment are excluded from the statistics. Dashed lines represent the eccentricities of 4 mm (green), 2 mm (red) and 1 mm (blue) distances to the brain-CSF boundary. Note the different scaling of the y-axes (Top row).

FIGURE 8

Figure 8. Accuracy comparison between CG- and DG-FEM for solving the MEG forward problem, i.e., the full B-field B (Equation 14), for different mesh resolutions. Visualized are the means (Left column) and the boxplots (Right column) of the RDM% (Top row) and MAG% (Bottom row), for tangentially oriented sources at logarithmically-scaled eccentricities. Dipoles not belonging to the brain compartment are excluded from the statistics. Dashed lines represent the eccentricities of 4 mm (green), 2 mm (red) and 1 mm (blue) distances to the brain-CSF boundary. Note the different scaling of the y-axes (Top row).

In this analysis the focus is on DG-FEM and the necessity of embedding the conservative flux (Equation 38) in the evaluation of the secondary B-field B^s. We will thus validate and compare the DG-FEM MEG forward methods with non-conservative (Equation 38) and conservative (Equation 39) flux.

The RDM% and MAG% statistical errors can be seen in Figure 5. Note that only a 4 mm mesh (seg_4_res_4) has been used, because the lower the resolution, the higher we expect the difference to be. With increasing source eccentricity, an overall increase of the RDM% (top row) and MAG% (bottom row) errors can be observed as shown by the mean error (left column) and by the boxplot statistics (right column). The boxplots indicate mainly increasing error statistics with regard to median, total range (TR), interquartile range (IQR) and also maxima for both conservative and non-conservative flux implementations. As a general result, the employment of the conservative flux (in green) delivers better results than the one of the non-conservative flux (in dark red). The difference between the two implementations is more evident with increasing eccentricity of the sources.

Let us now discuss in more detail the eccentricity of 0.9796, i.e., 1.59 mm from the brain-CSF boundary. Higher eccentricities do not have practical importance, as already explained in section 3.3. For the eccentricity of 0.9796, the maximum difference of 20 percentage points (pp) in mean RDM% is achieved between the conservative and the non-conservative DG flux approaches. For the least eccentric sources, this difference goes down to about 2 pp (see the 0.01 eccentricity in top left subfigure of Figure 5).

With regard to the boxplot of the RDM%, the median values of the conservative flux case are overall smaller than the ones of the non-conservative flux. For sources with eccentricity value of 0.9796 the RDM% median difference is greater than 20 pp; the IQR difference is approximately 15 pp and the TR is constant and similar for both approaches.

In the MAG% boxplot (right column), the much better performance of the conservative flux approach is especially clearly visible. The MAG% median difference reaches 40 pp for realistic sources of eccentricity 0.9796. For the same sources, the TRs, IQRs and means are in general large, with a ratio 1:4 between conservative and non-conservative flux values. For lower eccentricities, we observe overall smaller errors.

4.1.3. FEM Study 2: Convergence of DG Approach

Since we have seen in the last study that the conservative flux DG-FEM approach (Equation 39) performs remarkably better than the non-conservative approach (Equation 38), for the remainder of the paper, we proceed with DG-FEM as in Equation (49). The third study proposed is about the convergence of the DG-FEM for computing the secondary B-field B^s, when the mesh resolution is increased, namely from the coarsest resolution of 4 mm over 2 mm to the highest resolution of 1 mm. We studied the behavior of the RDM% and MAG% errors for 8,000 tangentially oriented and randomly distributed dipoles at different eccentricities. Results can be seen in Figure 6.

The RDM% and MAG% error mean curves (Figure 6, left column) are overall increasing with increasing source eccentricity, as hypothesized by the theory of the subtraction approach (Wolters et al., 2007) and well-known already from EEG results (Drechsler et al., 2009). Most importantly, for increasing mesh resolution, error statistics improve considerably. For the most relevant eccentricity of 0.9796, the highest resolved model (seg_1_res_1) reaches mean RDM% and MAG% errors of 1.5% and 0.1%, respectively. On the right column, we can study the boxplots of the RDM% and MAG% of the same scenario analyzed before. Both in the RDM% and MAG% cases, there is an overall increase of the median, TR and IQR when increasing the source eccentricity and decreasing the mesh resolution. If we focus on the 1 mm mesh and 0.9796 eccentricity, the RDM% median is only around 1.2%; the IQR is 0.8% and the TR reaches 20%. In particular, the IQR for dipoles of eccentricity 0.9796 increases drastically from 0.8% (1 mm) to almost 10% (2 mm) and 30% (4 mm). The TR behaves similarly. The median MAG% is extremely low, i.e., ≈0.017%; the IQR is ≈0.8% and the TR is ≈25%. For this eccentricity, we notice a huge difference among the three mesh resolutions: the medians grow from 0.017% (1 mm), to 2.8% (2 mm), up to 28.6% (4 mm). The same trend is noticeable for the IQR: 0.8% (1 mm), 7% (2 mm) and 60% (4 mm). However, these values are out of the displayed graph range. The TR again behaves similarly.

4.1.4. FEM Study 3: Comparison between CG and DG

The fourth analysis performed in this work is a comparison between CG- and DG-FEM for the MEG forward computation. RDM% and MAG%s are evaluated both for the secondary B-field B^s (Figure 7) and the full B-field B (Figure 8), following Equations (28) and (14), respectively.

In our following result discussion, we focus on the comparison between the two methods, rather than the performance of each method alone, which has been done for DG-FEM in section 4.1.3.

With regard to the secondary B-field B^s results, we will first analyze the mean RDM% curve (Figure 7, top left). In this plot we can distinguish the three different couples of curves: CG- and DG-FEM for 1 mm (seg_1_res_1), 2 mm (seg_2_res_2) and 4 mm (seg_4_res_4).

If we focus on the 1 mm analysis, we notice a high accuracy (up to around 1.5%) for eccentricities smaller or equal to 0.9796 (i.e., 1.59 mm from the CSF compartment). Even if in our current implementation, CG-FEM achieves slightly better results, the differences to DG-FEM are below 0.5 pp, so that in summary, DG-FEM constitutes an interesting alternative to the CG-FEM approach. Also for lower mesh resolutions of 2 and 4 mm, the performance of CG- and DG-FEM are very comparable for the realistic eccentricities up to 0.9796. A similar observation can be made for the mean MAG% curve, as the general trend for the three couples of curves (i.e., CG-DG 1 mm, CG-DG 2 mm, CG-DG 4 mm) is the same as before. When focusing on sources with eccentricity value of 0.9796, the mean MAG% difference between CG- and DG-FEM remains below 0.11 pp.

As for the boxplots for 1 mm mesh resolution (seg_1_res_1) and source eccentricity of 0.9796, the median RDM% difference is ≈0.4 pp (≈0.8 and ≈1.2% for CG- and DG-FEM, respectively); the IQR difference is around 0.2 pp (≈0.6 and ≈0.8% for CG- and DG-FEM, respectively) and the TR difference reaches almost 20 pp (Figure 7, top right). In the same scenario, the MAG% medians are identically extremely low, i.e., ≈0.015%. The IQRs also do not differ, while, again the TR difference is around 20 pp (Figure 7, bottom right).

The results when focusing on the full B-field B in Figure 8 are similar to the ones in Figure 7. Even for the full B-field, both the CG- and DG-FEM show an overall very high accuracy and a negligible difference, especially when focusing on the 1 mm study and source eccentricity of 0.9796. The mean RDM% (Figure 8, top left) is ≈0.9% for CG-FEM and ≈1.5% for DG-FEM; the mean MAG% is ≈ −0.015% for CG-FEM and ≈0.1% for DG-FEM (Figure 8, bottom right). With regard to the RDM% boxplot (Figure 8, top right), the medians are ≈0.8 and ≈1.15% for CG- and DG-FEM, respectively; the IQRs are ≈1 and ≈1.2% for CG- and DG-FEM, respectively, and the TRs are ≈3 and ≈20% for CG- and DG-FEM. In the MAG% boxplot (Figure 8, bottom right), we observe identical and extremely low values for the median (≈0.01%) and for the IQRs (≈0.8%). The difference of TRs is again bigger (≈20 pp) because of few outliers. Note that in Figure 8, we have omitted the errors for the lowest eccentricity of 0.01 because radial sources do not produce any magnetic field.

4.1.5. FEM Study 4: CG and DG in a Leaky Sphere Model

Motivated by the EEG results of Engwer et al. (2017), where DG-FEM could clearly outperform CG-FEM in skull leakage scenarios, this section is concerned with the comparison of CG- and DG-FEM for the same scenario, but for the MEG case. Therefore, a leaky sphere model (seg_2_res_2_r82) has been constructed using an outer skull radius of 82 mm (instead of 86 mm as in the previous sections), resulting in an only 2 mm thick spherical skull compartment. Then, a 2 mm resolution hexahedral model has been constructed, resulting in 10,080 skull leakages. Again, only tangentially-oriented dipoles have been examined.

In Figure 9, we computed RDM% and MAG% mean curves (left column) and boxplots (right column) for the leaky skull spherical model scenario (seg_2_res_2_r82) compared to the non-leaky skull spherical model scenario (seg_2_res_2).

FIGURE 9

Figure 9. Accuracy comparison for secondary B-field B^s computation (Equation 28) between CG-FEM (in warm colors) and DG-FEM with the conservative flux (in cold colors), in two different 2 mm hexahedral sphere models: seg_2_res_2 and seg_2_res_2_r82, described in Table 4. Visualized are the means (Left column) and the boxplots (Right column) of the RDM% (Top row) and MAG% (Bottom row), for tangentially oriented sources at logarithmically-scaled eccentricities. Dipoles not belonging to the brain compartment are excluded from the statistics. The dashed red line represents the eccentricity of 2 mm distance to the brain-CSF boundary. Note the different scaling of the y-axes (Top row).

We observe that, in contrast to the improvement that DG-FEM could achieve in the EEG case (Engwer et al., 2017), the skull leakages do not visibly influence the numerical simulation of the secondary B-field B^s and, since the primary B-field B^p is also not influenced, thereby also the full B-field and thus the MEG forward problem.

If we observe the plots in the left columns, we notice that the curves of the leaky scenarios are completely overlaying the curves of the non-leaky scenarios, both for CG- and DG-FEM and both for RDM% and MAG% mean curves.

Also in the boxplots we cannot distinguish the behavior of the RDM% and MAG% in the leaky or non-leaky scenarios.

4.2. DG-FEM MEG Study in a Realistic Head Model

In the last study, as a proof of concept, we simulated an auditory N1 MEG signal using the new DG-FEM method with conservative flux (Equation 39) in the 6 compartment realistically-shaped head volume conductor model. Following experimental evidence (Okamoto et al., 2007), the N1 current dipole was positioned in the secondary auditory cortex and oriented inwards-pointing and normally to the gray matter surface. The result is shown in Figure 10. The subfigure on the left represents a sagittal slice through the head model, color-coding the 6 tissue compartment with different conductivities. In the middle and right subfigures, the results for EEG and MEG forward problem are presented. More precisely, the dipolar electrical potential map with frontal negativity and right occipital positivity is visualized on a cropped volume of the hexahedral mesh together with the underlying source (black arrow). The normally-oriented B-field MEG results at the 275 magnetometers was interpolated and visualized, showing a dipolar pattern that is 90° rotated to the EEG one and, following the right-hand rule, the negativity (blue) is over central and the positivity over temporal areas, in line with the experimental results (Hämäläinen et al., 1993; Okamoto et al., 2007).

FIGURE 10

Figure 10. Exemplary EEG and MEG forward computation for an auditory source computed using DG-FEM in a realistically shaped head model. Hexahedral mesh with 2 mm resolution, 6 compartments, sagittal slice (Left); electric potential distribution visualized on the clipped volume conductor model in the sagittal plane where the auditory dipole (black cone) lies (Middle); MEG solution interpolated on the radial magnetometers including a volume rendering of the head model (Right).

5. Discussion

In this paper, we developed, implemented and evaluated one CG-FEM and two new DG-FEM approaches, a conservative and a non-conservative one, to solve the MEG forward problem. In section 2, we provided the mathematical theory for the CG-FEM and for the two new DG-FEM approaches with conservative and non-conservative discrete representation of the electrical flux. We started from the EEG formulation and continued with the MEG approaches. In section 3, we first described the implementation of the FEM-based MEG forward approaches in DUNEuro, a new modular C++ toolbox dedicated to solve partial differential equations in neuroscience. Furthermore, we presented the validation and evaluation platform for the new methods. In section 4, we presented the results of our analysis.

First, we tested and visualized the symmetry properties of the MEG analytical solution for a multi-layer homogeneous sphere model, as described in Remark 1 and as proven by Sarvas (1987). First of all, radial sources have a zero magnetic field outside a multi-layer sphere volume conductor model. Then, for tangential sources, only the primary B-field contributes to the full B-field for radial point-magnetometers (Figure 3A). For tangential sources and tangential point-magnetometers, we additionally showed that the more eccentric the source is, the more its primary B-field contributes to the full B-field relative to the contribution of the secondary B-field. The deeper the tangential source is, the more the secondary B-field weakens the primary B-field until for sources in the middle of the sphere model, where the primary and secondary B-fields totally compensate for each other (Figure 4 together with Figure 3B). In contrast, the more eccentric the sources are, the less symmetric the return currents are and the less their secondary B-field compensates the magnetic field of the primary current (Figure 4). This is in line with the fact that the strength of the full B-field for both radial and tangential magnetometers is decreasing for decreasing eccentricities (Figure 3B), as expected by the theory where it is proven that MEG sensors are blind for radial sources (Sarvas, 1987).

In a second analysis, we studied how large the influence of a conservative representation of the electrical flux in the computation of the secondary B-field is by adopting the DG-FEM is adopted. By comparing the DG-FEM with a conservative (Equation 39) and non-conservative (Equation 38) flux in a 4 mm multi-layer homogeneous sphere model, the high importance of DG-FEM with conservative flux could be worked out, outperforming the non-conservative DG-FEM scheme in all cases (Figure 5). In light of these results, the conservative flux DG-FEM was then used in all consecutive studies.

Results from the third study show the convergence of the DG-FEM numerical solutions toward the analytical solution when the resolution of the meshes is increased from 4 mm over 2 mm down to 1 mm (Figure 6).

From our comparison studies between CG- and DG-FEM regarding the secondary and the full B-fields, we first of all noticed that the accuracy for the 1 mm mesh resolution is extremely accurate for both methods: the mean RDM% is only up to ≈1.5% and the mean MAG% only up to ≈0.1% for sources with realistic eccentricities of 0.9796 (i.e., 1.59 mm to the next conductivity jump at the brain-CSF boundary) (Figures 7, 8). With ≈1.15% median RDM% and ≈0.01% MAG% for DG-FEM for the 1 mm mesh, also the IQRs are very low. The result is only slightly influenced by some few outliers, which showed higher numerical errors. However, these errors might be avoided by better controlling the source position with regard to the mesh, for example by only allowing sources in the center of the cortical hexahedra.

To the best of our knowledge, not many recent studies on finite element methods applied to solve the MEG forward problem have been presented, and none of them about DG-FEM. (Van den Broek et al., 1996) applied a CG-FEM approach in a 10 cm single-layer homogeneous sphere model and RDM% errors were measured. The minimum RDM% found by these authors for far less eccentric sources of 0.95 was 3%. Still, the comparison is not straightforward because of the different approaches, the different element meshes (tetrahedrons vs. hexahedrons), and the different source models which have been used. In general, tetrahedral meshes can better approximate surfaces but, for realistic head models, the generation of such models is difficult in practice and might cause unrealistic model features, e.g., holes in tissue compartments such as the foramen magnum and the optic canals in the skull are often artificially closed to allow constrained Delaunay tetrahedralization (CDT). Furthermore, CDT modeling necessitates the generation of nested, non-intersecting, and non-touching surfaces. However, in reality, surfaces might touch, for example, the inner skull and the outer brain surfaces. Hexahedral models, as investigated here, have larger geometry approximation errors, but do not suffer from the above limitations and can be easily generated from voxel-based MRI data. However, with new methods like in Nüßing et al. (2016), such geometry approximation errors can be avoided without the need of generating geometry conforming tetrahedral meshes.

van den Broek et al. (1998) used both a sphere model and a realistically shaped model. In both cases only three compartments were modeled, namely the brain, the skull and the scalp (brain and scalp with the same conductivity values, and skull with a 1:80 conductivity ratio). The authors used a lower amount of sources and lower mesh resolutions, but locally-refined tetrahedra meshes. In both scenarios, magnetometer sensors were covering only the top half of the models.

A CG-FEM MEG forward modeling study in a human (and rabbit) head volume conductor model was performed by Haueisen et al. (1995). The authors distinguished 12 or more homogeneous and isotropic realistically shaped head tissue compartments and used 2 mm FEM models. Since the focus was on sensitivity analysis and suppression ratio (i.e., the magnetic field of radial dipole divided by the one of the corresponding tangential dipole, was found to be in average 0.19 ± 0.07 in the realistic human head model) and not on validation in sphere models like in our study, we can not further compare these results to our results.

Another example of a CG-FEM and Biot-Savart's law scheme used to compute the electric potential and the B-field was presented by Schimpf et al. (2002). Similar to our approach, the authors used a 1 mm hexahedral mesh of a 4-layer piecewise homogeneous and isotropic sphere model. Also the arrangement of sources and sensors was similar to our work. The main focus of their work was, however, on source modeling: it was found that from the different tested source modeling approaches, the subtraction approach, also used in our study at hand, was the most accurate one.

In Vorwerk (2016, Ch. 2.10.4), for solving the MEG forward problem, three different CG-FEM source modeling approaches (i.e., subtraction, Saint Venant and partial integration) were compared in a 1 mm hexahedral (and in tetrahedral) meshes. Both the secondary and the full B-fields were examined against the analytical solution in a multi-layer homogeneous sphere model for tangentially oriented magnetometers. Also in this comparison it was found that the subtraction approach outperforms the other source modeling methods with regard to numerical accuracy for all sources apart from the most eccentric one. The subtraction method is therefore most sensitive to very close conductivity jumps and thus needs high resolution meshes especially in the source area, a result which is in line with ours. Deeper comparisons are again not easy because of different set-ups, but we are planning a direct comparison of the SimBio³ code used by Vorwerk (2016) and our DUNEuro implementation in future studies.

In Vorwerk et al. (2014), a guideline for EEG and MEG forward modeling using CG-FEM Saint Venant modeling was presented in realistic head models with a varying number of layers and conductivity profiles. The main result was that it is highly recommended to include the CSF and distinguish between gray and white matter and that, especially for the MEG, the modeling of skull spongiosa and compacta might be neglected. Furthermore, the numerical errors of a lower resolved (about 1 million nodes) 6 compartment anisotropic (6CA) model in reference to a higher resolved (about 2 millions nodes) version of 6CA were studied and expressed in terms of topography and magnitude errors: 95% of the sources had an RDM% of less than 2.5% and a MAG% of less than 10%.

Our last study was about the influence of leaky points on the computation of the secondary B-field when DG-FEM is adopted (Figure 9). In this analysis we considered two different multi-layer homogeneous sphere models, namely seg_2_res_2 and seg_2_res_2_r82. The difference between the two models is that in seg_2_res_2_r82 the thickness of the skull compartment is deliberately reduced so that 10,080 leaky points are present. When comparing CG- and DG-FEM in the leaky model (i.e., seg_2_res_2_r82) and in the non-leaky model (i.e., seg_2_res_2), we observed that the results of the computation of the secondary B-fields are almost identical. This means that the skull leakages neither cause additional MEG forward modeling errors for DG-FEM, nor for CG-FEM. The situation is thus different from the EEG case, where remarkable errors for CG-FEM forward modeling were shown, while DG-FEM could strongly alleviate these additional leakage errors (Engwer et al., 2017). For MEG, in case of tangential sources, the return currents mainly flow parallel to the inner skull surface in the close environment of the source, so that the leakages do not affect the overall MEG forward solution. We have to underline the fact that the results obtained in leaky scenarios are not to be confused with those where real holes of a certain diameter, e.g., from trepanation, are present in the skull compartment. The skull leakages investigated in Engwer et al. (2017) and, consequently, in this work are due to erroneous or, in general, poor representation of the skull compartment and not to real holes in the skull compartment. Lau et al. (2014) found that MEG signals are influenced by skull defects such as post-surgical skull openings. They examined the influence of skull holes in MEG signals via in vivo rabbit brains experiments, finding that the MEG signal amplitude reduced by as much as 20%, especially if the source is central under the skull defect. Their conclusion is that MEG source modeling requires realistic volume conductor head models that incorporate skull defects. Furthermore, Lau et al. (2016) showed that also MEG inverse solutions are affected by skull defects. In particular, ignoring skull defects in the head model during reconstruction displaced and reoriented sources under a skull defect, and when skull defects were incorporated in the head model with their physical conductivity, the location and orientation errors were mostly eliminated.

A further important aspect to discuss is that, if a combined EEG and MEG source reconstruction is strived for (Fuchs et al., 1998; Aydin et al., 2015), the same forward model should be used for both EEG and MEG, because of considerable advantages in terms of implementation, accuracy and computational cost efficiency, as the MEG forward model is also based on the electric potential and thus the numerical solution of the EEG forward problem. We therefore employed the same method (CG- or DG-FEM, with conservative or non-conservative flux representation) for both EEG and MEG in our work at hand. Accordingly, in case of EEG or combined MEG/EEG source reconstruction in possibly leaky head models (e.g., in temporal bone areas or, more generally, in children investigations), the usage of DG-FEM is recommended. In fact, DG-FEM clearly improves EEG forward solutions in leaky models (Engwer et al., 2017) and, at the same time, delivers reliable and accurate MEG solutions, as shown in the study at hand.

In this study, we did not evaluate the computational costs of the CG- and DG-FEM schemes for the computation of the MEG forward solution. Because of the higher number of degrees of freedom, DG-FEM is computationally more expensive than CG-FEM. However, the FEM transfer matrix approach (section 2.4) considerably reduces the computational costs of both approaches, so that this aspect gets less relevant for practical applications.

We now discuss possibilities for further accuracy increase that we plan to evaluate in our future work. In this study, sources were just chosen randomly, i.e., the influence of the source position relative to an element of the discretization was not yet investigated. It is well known that the combination of computing leadfields only for the most accurate sources combined with leadfield inter- and extrapolation techniques for other sources might not only speed up computations, but might also further increase numerical accuracy (Yvert et al., 2001; Vorwerk, 2011). In the DG-FEM scheme, indeed, already in the EEG forward computation (see Equation 26), the contribution given by the integral over Γ_int can reach high values when the source is relatively close to a quadrature point on the internal skeleton, because of the singularity in ∇u^∞. Moreover, in Drechsler et al. (2009), an analytical expression for ∇u^∞ was derived for isotropic and anisotropic conductivity distributions in the source space. A further future goal will thus be its implementation and use to further decrease the numerical errors in our FEM implementations, both on the CG- and DG-FEM sides. In addition, the degrees of polynomials in $V_h^1$ can be increased, together with the order of the Raviart-Thomas function space used to extend the conservative flux into the volume of each element. On the other hand, increasing the order of function spaces results in increased computational costs, so this intervention should be treated carefully. Furthermore, the DG-FEM constitutes the first step for the UDG-FEM implementation. This method, already tested in an EEG study (Nüßing et al., 2016), reduces the geometrical error of the forward simulations in hexahedral models while drastically decreasing the computational cost and thus its application to the MEG forward modeling represents an interesting future goal.

Overall the newly implemented conservative flux DG-FEM scheme offers an interesting new EEG and MEG forward modeling approach. It can be used especially in leakage scenarios and, in general, for comparison purposes, not only in EEG and MEG source analysis, but also in bioelectromagnetism applications, i.e., including also the simulation of transcranial electric and/or magnetic stimulation (Miranda et al., 2006; Datta et al., 2013; Windhoff et al., 2013; Wagner et al., 2014).

6. Conclusions

We presented theory, validation and evaluation of three finite element method (FEM) approaches for the MEG forward problem, namely the continuous Galerkin FEM (CG-FEM), as well as two new approaches, the discontinuous Galerkin FEM (DG-FEM) with a conservative and a non-conservative flux implementation. All three methods have been implemented in the DUNEuro software module. Statistical validations and evaluations have been performed on multi-layer homogeneous sphere models represented via hexahedral meshes and the subtraction approach has been adopted as source model. DG-FEM with conservative flux implementation, i.e., a main feature of a DG-FEM discretization, turned out to be superior to the non-conservative flux variant. The new DG-FEM method showed proper convergence behavior with increasing mesh resolution. When compared to the CG-FEM, DG-FEM provided results that are in a comparable range of high accuracy. Furthermore, both methods are able to model realistic head volume conductor models with their tissue inhomogeneities and anisotropies. In contrast to EEG studies, the so-called skull leakage effects did not play a crucial role for MEG. However, for EEG or combined MEG/EEG source analysis scenarios, DG-FEM offers an interesting new alternative to CG-FEM, considering the importance of a high accuracy of the forward problem solution in MEG/EEG source reconstruction. Finally, the DG-FEM MEG forward simulation in a realistic head model for an auditory source resulted in EEG and MEG topographies that are in line with practical findings in the field of auditory evoked responses.

Author Contributions

MP, CE, and CW conceived the study. MP wrote the code with the supervision of AN and CE. JV and CW provided spherical and realistic head models. MP and JV constructed the FEM models. MP performed the simulations and the analysis, produced the images, interpreted the results, and wrote the paper. All authors took part in the scientific discussion at multiple stages of the study and provided feedback from the modeling (CW, CE, HB, RO), theoretical (CE, CW, AN, JV), and technical (AN, CE) perspective. All authors reviewed the manuscript and approved it for publication.

Funding

This work was supported by EU project ChildBrain (Marie Curie Innovative Training Networks, grant agreement 641652), by the Deutsche Forschungsgemeinschaft (DFG), project WO1425/7-1, by the Priority Program 1665 of the DFG (WO1425/5-2), and partially by the Cluster of Excellence 1003 of the Deutsche Forschungsgemeinschaft (DFG EXC 1003 Cells in Motion).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We wish to thank Professor Marco M. Fato and Dr. Gabriele Arnufo (University of Genova) for their support.

Footnotes

1. ^http://www.dune-project.org

2. ^www.duneuro.org

3. ^https://www.mrt.uni-jena.de/simbio/index.php/Main_Page

References

Acar, Z. A., and Makeig, S. (2010). Neuroelectromagnetic forward head modeling toolbox. J. Neurosci. Methods 190, 258–270. doi: 10.1016/j.jneumeth.2010.04.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Alkämper, M., Dedner, A., Klöfkorn, R., and Nolte, M. (2016). The DUNE-ALUGrid module. Arch. Numer. Softw. 4, 1–28.

Awada, K. A., Jackson, D. R., Williams, J. T., Wilton, D. R., Baumann, S. B., and Papanicolaou, A. C. (1997). Computational aspects of finite element modeling in EEG source localization. IEEE Trans. Biomed. Eng. 44, 736–752. doi: 10.1109/10.605431

PubMed Abstract | CrossRef Full Text | Google Scholar

Aydin, Ü., Rampp, S., Wollbrink, A., Kugel, H., Cho, J.-H., Knösche, T. R., et al. (2017). Zoomed MRI guided by combined EEG/MEG source analysis: a multimodal approach for optimizing presurgical epilepsy work-up and its application in a multi-focal epilepsy patient case study. Brain Topogr. 30, 417–433. doi: 10.1007/s10548-017-0568-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Aydin, Ü., Vorwerk, J., Dümpelmann, M., Küpper, P., Kugel, H., Heers, M., et al. (2015). Combined EEG/MEG can outperform single modality EEG or MEG source reconstruction in presurgical epilepsy diagnosis. PLoS ONE 10:e0118753. doi: 10.1371/journal.pone.0118753

PubMed Abstract | CrossRef Full Text | Google Scholar

Bastian, P., Blatt, M., Dedner, A., Engwer, C., Klöfkorn, R., Kornhuber, R., et al. (2008a). A generic grid interface for parallel and adaptive scientific computing. Part II: implementation and tests in DUNE. Computing 82, 121–138. doi: 10.1007/s00607-008-0004-9

CrossRef Full Text | Google Scholar

Bastian, P., Blatt, M., Dedner, A., Engwer, C., Klöfkorn, R., Ohlberger, M., et al. (2008b). A generic grid interface for parallel and adaptive scientific computing. Part I: abstract framework. Computing 82, 103–119. doi: 10.1007/s00607-008-0003-x

CrossRef Full Text | Google Scholar

Bastian, P., Heimann, F., and Marnach, S. (2010). Generic implementation of finite element methods in the distributed and unified numerics environment (DUNE). Kybernetika 46, 294–315.

Google Scholar

Baumann, S. B., Wozny, D. R., Kelly, S. K., and Meno, F. M. (1997). The electrical conductivity of human cerebrospinal fluid at body temperature. IEEE Trans. Biomed. Eng. 44, 220–223. doi: 10.1109/10.554770

PubMed Abstract | CrossRef Full Text | Google Scholar

Bertrand, O., Thevenet, M., and Perrin, F. (1991). “3-D finite element method in brain electrical activity studies,” in Biomagnetic Localization and 3D Modeling, Report TKK-F-A689, eds J. Nenonen, H.-M. Rajala, and T. Katila (Espoo: Helsinki University of Technology, Department of Technical Physics, Laboratory of Biomedical Engineering).

Blatt, M., and Bastian, P. (2007). “The iterative solver template library,” in Applied Parallel Computing. State of the Art in Scientific Computing, Lecture Notes in Computer Science, eds B. Kågström, E. Elmroth, J. Dongarra, and J. Waśniewski (Berlin; Heidelberg: Springer), 666–675.

Google Scholar

Brette, R., and Destexhe, A. (2012). Handbook of Neural Activity Measurement. New York, NY: Cambridge University Press. doi: 10.1017/CBO9780511979958

CrossRef Full Text | Google Scholar

Cook, M. J., and Koles, Z. J. (2006). “A high-resolution anisotropic finite-volume head model for EEG source analysis,” in Engineering in Medicine and Biology Society, 2006. EMBS'06. 28th Annual International Conference of the IEEE (New York, NY: IEEE), 4536–4539.

PubMed Abstract | Google Scholar

Dannhauer, M., Lanfer, B., Wolters, C. H., and Knösche, T. R. (2011). Modeling of the human skull in EEG source analysis. Hum. Brain Mapp. 32, 1383–1399. doi: 10.1002/hbm.21114

PubMed Abstract | CrossRef Full Text | Google Scholar

Datta, A., Zhou, X., Su, Y., Parra, L. C., and Bikson, M. (2013). Validation of finite element model of transcranial electrical stimulation using scalp potentials: implications for clinical dose. J. Neural Eng. 10:036018. doi: 10.1088/1741-2560/10/3/036018

PubMed Abstract | CrossRef Full Text | Google Scholar

De Munck, J. C., Van Dijk, B. W., and Spekreijse, H. (1988). Mathematical dipoles are adequate to describe realistic generators of human brain activity. IEEE Trans. Biomed. Eng. 35, 960–966. doi: 10.1109/10.8677

PubMed Abstract | CrossRef Full Text | Google Scholar

Drechsler, F., Wolters, C. H., Dierkes, T., Si, H., and Grasedyck, L. (2009). A full subtraction approach for finite element method based source analysis using constrained Delaunay tetrahedralisation. NeuroImage 46, 1055–1065. doi: 10.1016/j.neuroimage.2009.02.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Engwer, C., Vorwerk, J., Ludewig, J., and Wolters, C. H. (2017). A discontinuous galerkin method to solve the EEG forward problem using the subtraction approach. SIAM J. Sci. Comput. 39, B138–B164. doi: 10.1137/15M1048392

CrossRef Full Text | Google Scholar

Fortin, M., and Brezzi, F. (1991). Mixed and Hybrid Finite Element Methods. New York, NY: Springer-Verlag.

Google Scholar

Fuchs, M., Wagner, M., Wischmann, H.-A., Köhler, T., Theißen, A., Drenckhahn, R., et al. (1998). Improving source reconstructions by combining bioelectric and biomagnetic data. Electroencephalogr. Clin. Neurophysiol. 107, 93–111. doi: 10.1016/S0013-4694(98)00046-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Gramfort, A., Papadopoulo, T., Olivi, E., and Clerc, M. (2011). Forward field computation with OpenMEEG. Comput. Intell. Neurosci. 2011:923703. doi: 10.1155/2011/923703

PubMed Abstract | CrossRef Full Text | Google Scholar

Hämäläinen, M., Hari, R., Ilmoniemi, R. J., Knuutila, J., and Lounasmaa, O. V. (1993). Magnetoencephalography—theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Mod. Phys. 65:413. doi: 10.1103/RevModPhys.65.413

CrossRef Full Text | Google Scholar

Haueisen, J., Ramon, C., Czapski, P., and Eiselt, M. (1995). On the influence of volume currents and extended sources on neuromagnetic fields: a simulation study. Ann. Biomed. Eng. 23, 728–739. doi: 10.1007/BF02584472

PubMed Abstract | CrossRef Full Text | Google Scholar

Ilmoniemi, R. (1995). Radial anisotropy added to a spherically symmetric conductor does not affect the external magnetic field due to internal sources. Europhys. Lett. 30:313. doi: 10.1209/0295-5075/30/5/011

CrossRef Full Text | Google Scholar

Lau, S., Flemming, L., and Haueisen, J. (2014). Magnetoencephalography signals are influenced by skull defects. Clin. Neurophysiol. 125, 1653–1662. doi: 10.1016/j.clinph.2013.12.099

PubMed Abstract | CrossRef Full Text | Google Scholar

Lau, S., Güllmar, D., Flemming, L., Grayden, D. B., Cook, M. J., Wolters, C. H., et al. (2016). Skull defects in finite element head models for source reconstruction from magnetoencephalography signals. Front. Neurosci. 10:141. doi: 10.3389/fnins.2016.00141

PubMed Abstract | CrossRef Full Text | Google Scholar

Marin, G., Guerin, C., Baillet, S., Garnero, L., and Meunier, G. (1998). Influence of skull anisotropy for the forward and inverse problem in EEG: simulation studies using FEM on realistic head models. Hum. Brain Mapp. 6, 250–269. doi: 10.1002/(SICI)1097-0193(1998)6:4<250::AID-HBM5>3.0.CO;2-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Meijs, J. W., Weier, O. W., Peters, M. J., and Van Oosterom, A. (1989). On the numerical accuracy of the boundary element method (EEG application). IEEE Trans. Biomed. Eng. 36, 1038–1049. doi: 10.1109/10.40805

PubMed Abstract | CrossRef Full Text | Google Scholar

Miranda, P. C., Lomarev, M., and Hallett, M. (2006). Modeling the current distribution during transcranial direct current stimulation. Clin. Neurophysiol. 117, 1623–1629. doi: 10.1016/j.clinph.2006.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Montes-Restrepo, V., van Mierlo, P., Strobbe, G., Staelens, S., Vandenberghe, S., and Hallez, H. (2014). Influence of skull modeling approaches on EEG source localization. Brain Topogr. 27, 95–111. doi: 10.1007/s10548-013-0313-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Mosher, J. C., Leahy, R. M., and Lewis, P. S. (1999). EEG and MEG: forward solutions for inverse methods. IEEE Trans. Biomed. Eng. 46, 245–259. doi: 10.1109/10.748978

PubMed Abstract | CrossRef Full Text | Google Scholar

Murakami, S., and Okada, Y. (2006). Contributions of principal neocortical neurons to magnetoencephalography and electroencephalography signals. J. Physiol. 575, 925–936. doi: 10.1113/jphysiol.2006.105379

PubMed Abstract | CrossRef Full Text | Google Scholar

Nédélec, J.-C. (1980). Mixed finite elements in r³. Numerische Mathematik 35, 315–341. doi: 10.1007/BF01396415

CrossRef Full Text | Google Scholar

Nüßing, A., Wolters, C. H., Brinck, H., and Engwer, C. (2016). The unfitted discontinuous Galerkin method for solving the EEG forward problem. IEEE Trans. Biomed. Eng. 63, 2564–2575. doi: 10.1109/TBME.2016.2590740

PubMed Abstract | CrossRef Full Text | Google Scholar

Okamoto, H., Stracke, H., Wolters, C. H., Schmael, F., and Pantev, C. (2007). Attention improves population-level frequency tuning in human auditory cortex. J. Neurosci. 27, 10383–10390. doi: 10.1523/JNEUROSCI.2963-07.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Pursiainen, S., Vorwerk, J., and Wolters, C. H. (2016). Electroencephalography (EEG) forward modeling via H (div) finite element sources with focal interpolation. Phys. Med. Biol. 61:8502. doi: 10.1088/0031-9155/61/24/8502

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramon, C., Schimpf, P., Haueisen, J., Holmes, M., and Ishimaru, A. (2004). Role of soft bone, CSF and gray matter in EEG simulations. Brain Topogr. 16, 245–248. doi: 10.1023/B:BRAT.0000032859.68959.76

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarvas, J. (1987). Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. Phys. Med. Biol. 32:11. doi: 10.1088/0031-9155/32/1/004

PubMed Abstract | CrossRef Full Text | Google Scholar

Schimpf, P. H., Ramon, C., and Haueisen, J. (2002). Dipole models for the EEG and MEG. IEEE Trans. Biomed. Eng. 49, 409–418. doi: 10.1109/10.995679

PubMed Abstract | CrossRef Full Text | Google Scholar

Stenroos, M., and Sarvas, J. (2012). Bioelectromagnetic forward problem: isolated source approach revis (it) ed. Phys. Med. Biol. 57:3517. doi: 10.1088/0031-9155/57/11/3517

PubMed Abstract | CrossRef Full Text | Google Scholar

Van den Broek, S., Zhou, H., and Peters, M. (1996). Computation of neuromagnetic fields using finite-element method and Biot-Savart law. Med. Biol. Eng. Comput. 34, 21–26. doi: 10.1007/BF02637018

PubMed Abstract | CrossRef Full Text | Google Scholar

van den Broek, S. P., Reinders, F., Donderwinkel, M., and Peters, M. (1998). Volume conduction effects in EEG and MEG. Electroencephalogr. Clin. Neurophysiol. 106, 522–534. doi: 10.1016/S0013-4694(97)00147-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Vatta, F., Meneghini, F., Esposito, F., Mininel, S., and Di Salle, F. (2009). Solving the forward problem in EEG source analysis by spherical and fdm head modeling: a comparative analysis-biomed 2009. Biomed. Sci. Instrum. 45, 382–388.

PubMed Abstract | Google Scholar

Vorwerk, J. (2011). Comparison of Numerical Approaches to the EEG Forward Problem. Diploma Thesis in Mathematics, Westfälische Wilhelms-Universität Münster.

Vorwerk, J. (2016). New Finite Element Methods to Solve the EEG/MEG Forward Problem. Ph.D. thesis, Ph.D. thesis in Mathematics, Westfälische Wilhelms-Universität Münster.

Vorwerk, J., Cho, J.-H., Rampp, S., Hamer, H., Knösche, T. R., and Wolters, C. H. (2014). A guideline for head volume conductor modeling in EEG and MEG. NeuroImage 100, 590–607. doi: 10.1016/j.neuroimage.2014.06.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, S., Lucka, F., Vorwerk, J., Herrmann, C. S., Nolte, G., Burger, M., et al. (2016). Using reciprocity for relating the simulation of transcranial current stimulation to the EEG forward problem. Neuroimage 140, 163–173. doi: 10.1016/j.neuroimage.2016.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, S., Rampersad, S., Aydin, Ü., Vorwerk, J., Oostendorp, T., Neuling, T., et al. (2014). Investigation of tDCS volume conduction effects in a highly realistic head model. J. Neural Eng. 11:016002. doi: 10.1088/1741-2560/11/1/016002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wendel, K., Narra, N. G., Hannula, M., Kauppinen, P., and Malmivuo, J. (2008). The influence of CSF on EEG sensitivity distributions of multilayered head models. IEEE Trans. Biomed. Eng. 55, 1454–1456. doi: 10.1109/TBME.2007.912427

PubMed Abstract | CrossRef Full Text | Google Scholar

Windhoff, M., Opitz, A., and Thielscher, A. (2013). Electric field calculations in brain stimulation based on finite elements: an optimized processing pipeline for the generation and usage of accurate individual head models. Hum. Brain Mapp. 34, 923–935. doi: 10.1002/hbm.21479

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolters, C. H., Grasedyck, L., and Hackbusch, W. (2004). Efficient computation of lead field bases and influence matrix for the FEM-based EEG and MEG inverse problem. Inverse Prob. 20:1099. doi: 10.1088/0266-5611/20/4/007

CrossRef Full Text | Google Scholar

Wolters, C. H., Köstler, H., Möller, C., Härdtlein, J., and Anwander, A. (2007). “Numerical approaches for dipole modeling in finite element method based source analysis,” in International Congress Series, vol. 1300 (Elsevier), 189–192.

Google Scholar

Yvert, B., Crouzeix-Cheylus, A., and Pernier, J. (2001). Fast realistic modeling in bioelectromagnetism using lead-field interpolation. Hum. Brain Mapp. 14, 48–63. doi: 10.1002/hbm.1041

PubMed Abstract | CrossRef Full Text | Google Scholar