On the Geometric Sensitivity of the MEG Inversion Algorithm

Magnetoencephalography (MEG) is a non-invasive technique that measures the magnetic fields produced by brain's electrical currents. The basic inverse problem of magnetoencephalography (MEG) consists in estimating the neuronal current in the brain from the measurement of the magnetic field outside the head. The relative inversion algorithms, existing in medical devices, for the identification of excitation sources inside the brain using MEG data, are based on the assumption that the geometry of the brain-head system is spherical. Here, we present the error of identification of the dipole (r0,Q)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${({\varvec{r}}}_{0}, {\varvec{Q}})$$\end{document} using the inversion algorithm of spherical geometry. In particular, taking magnetoencephalographic measurements from realistic ellipsoidal model and using these data in a spherical model lead to a structural error. For the purpose of this paper, we make numerical investigation for different positions of the source, especially, locating sources at different depths and see how the error depends on depth. These results may have useful implications for the interpretation of the reconstructions obtained via the existing approaches.


Introduction
The human brain is the most complex-organized structure known to exist, and, for us, it is also the most important. There are at least 10 10 neurons in the outermost layer of the brain and the cerebral cortex. The electromagnetic brain activity is generated via instantaneous flux of ions in the neurons as a result of intrinsic electrochemical reactions [1]. The neuronal current produces an electric potential which can be measured on the scalp, as well as a small magnetic flux. Measurements of the electric potential on the scalp lead to the Electroencephalography (EEG), while measurements of the magnetic flux outside the head lead to the magnetoencephalography (MEG), which are two of the most important techniques for studying brain neuronal activity in vivo. For both of these modalities, we know the direct mathematical problems of determining the electric potential and the magnetic flux when the interior neuronal sources, and the geometry of the brain-head system, are given, as well as the inverse mathematical problems of identifying the sources once the measured EEG and MEG data are obtained [2][3][4][5]. Hermann von Helmholtz showed in 1853 that this problem has no unique solution. During the last thirty years, direct and inverse problems of EEG and MEG have been studied extensively [1]. The geometry of the head affects the solutions of the mathematical problems of the electromagnetic brain activity since the ellipsoidal shape brings into the problem another more severe difficulty, since solving boundary value problems in ellipsoidal geometry is by far not a trivial analytic procedure. A triaxial ellipsoid with average semi-axes equal to 6 cm, 6.5 cm, and 9 cm provides an approximation of the average human brain which is much better than the most used spherical model. The first effort to solve the direct EEG problems appeared in [6,7]. Ever since a number of solutions for the EEG and MEG problems, with dipole as well as with continuous source distributions in ellipsoidal geometry, have been published [1,[8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24].
The machinery used in EEG and MEG practical diagnostics uses algorithms which are based on the assumption that the head has spherical shape [25,26]. Consequently, they record data from an ellipsoidal shape of brain and interpret them as they were coming from a sphere. Thus, an error occurs on the estimation of the location and the moment of a dipole. Since the solutions for the inverse MEG dipole problem in spherical and ellipsoidal geometry were known, the challenge was to transfer information from one geometrical system to the other, and anyone who has tried to do something like that immediately recognizes the difficulty of such a task [27]. We need to represent elements of a Hilbert space in spherical and ellipsoidal backgrounds although general representations are not possible at the analytic level [25]. However, as we demonstrate in this work, after some long and tedious calculations, these representations were achieved at least for the needed low degree eigen functions. The final result is the dependence of the error on the values of the principal eccentricities of the ellipsoid.
The paper is organized as follows. In Sect. 2, we provide a brief formulation of the mathematical problem, as well as the necessary information from the theory of ellipsoidal harmonics. That will make the paper readable without having to look at other complementary sources. Section 3 involves the general expansions of the electric potential in spherical and ellipsoidal harmonics, as well as the appropriate connection formulae. This is the core of the present work. For the sake of completeness, we provide a brief proof of Meg inversion algorithm [1,18,26,28] in Sects. 4 and 5. The main results are exposed, and Sect. 6 deals with the numerical implementation of the solution obtained results. The paper closes with concluding remarks in Sect. 7.

Mathematical Formulation
Let Ω be a bounded, connected, and homogeneous conductor with a smooth boundary S ∂Ω. The domain Ω provides a simplified geometrical model of the brain as an isotropic and homogeneous conductor with conductivity σ. The domain Ω is surrounded by the non-conductive exterior space c . The magnetic permeability μ 0 is assumed to be the same both in Ω and in c . A primary neuronal current J P with support in generates the electric and magnetic activity.
The mathematical formulations of the brain-imaging techniques of electroencephalography (EEG) and magnetoencephalography (MEG) are based on the quasistatic theory of electromagnetism. The quasi-static approximation of Maxwell's equations [1,26,[29][30][31][32][33][34][35][36] reads ∇ · E(r) 0, where E and B denote the electric and magnetic flux fields, respectively [27,37]. Since E is irrotational, it can be represented by an electric potential u, such that In the interior of a homogeneous conductor , a local neuronal excitation is represented by the equivalent dipole current: where δ stands for the Dirac measure at a fixed point r 0 with a dipole moment equal to Q. The primary current J p (r) induces an electric field E in the interior conductive space, which in turn generates an inductive volume current with density J v (r) σ • E(r) resulting to the total current density: The current J generates an electromagnetic wave, which propagates in the interior as well as in the exterior to the conducting space.
The potential u + is the field recorded in any electroencephalogram. In particular, we denote the electric potential in the interior space Ω, by u − and in the exterior space c , by u + . Combining Eqs. (2), (5), and (7), we obtain the Poisson equation: which the interior potential u-must satisfy in Ω.
In the source-free space c , the potential u + solves the Laplace equation: On the surface S, the following transmission conditions hold: where the outward normal differentiation on the surface is considered. Conditions (10) and (11) state the continuity of the potential function as well as the continuity of the normal component of current density on S.
In addition, we assume the asymptotic behavior at infinity: The relationship between dipole (r 0 , Q) and B is expressed by the Geselowitz equation:

Ellipsoidal System
A complete analysis of the material presented in this section can be found in [25]. Let S e denote the triaxial ellipsoid for which rectangular coordinates are provided as follows: where α 1 , α 2 , α 3 are three fixed parameters determining the reference semi-axes.
The basic ellipsoid (14) introduces an ellipsoidal system with coordinates (ρ, μ, ν) and semi-focal distances h 1 ,h 2 ,h 3 where The ellipsoidal coordinates (ρ, μ, ν) involve the ellipsoidal variable ρ ∈ [h 2 , +∞) and the hyperboloidal variables The connection between the ellipsoidal and the Cartesian systems is given by The coordinate ρ plays the role of the radial variable r, while μ and ν correspond to the angular variable θ and ϕ in spherical coordinates.
The ellipto-spherical coordinate system (ρ, θ ε , ϕ ε ) combines the ellipsoidal variable that specifies the family of confocal ellipsoids with the eccentric angular variables of the spherical system. It is defined by In ellipsoidal coordinates, the surface S e given in (14) corresponds to ρ α 1 , and it represents the boundary of the brain. The interior space e is defined by the interval ρ ∈ [h 2 , a 1 ), and it is characterized by the conductivity σ. The exterior non-conductive space c e is defined by ρ ∈ (a 1 , ∞). Separation of variables for Laplace's Eq. (9) in the ellipsoidal coordinate system leads to the Lamé equation: for each one of the factors E(ρ), E(μ), and E(ν) that form the interior harmonic function: In Eq. (18), the parameters P and n are constants that define, in a complicated way, the degree n and the order m of the interior ellipsoidal harmonic (19).
The products E m n (μ), E m n (ν) defined on the surface of any specific ellipsoid, are known as surface ellipsoidal harmonics, and they form a complete orthogonal set of surface eigen functions.
We define the normalization constants γ m n as follows: where the employed symbol of integration indicates integration over the surface of the ellipsoid ρ α 1 .
In the present work, the following normalization constants are to be used:

Spherical System
Next, we restrict our consideration to the case where the domain s is a sphere of radius α centered at the origin. The domain s is bounded by the spherical surface S s for which rectangular coordinates are provided as follows: In spherical coordinates, the surface S s given in (30) corresponds to r α and it represents the boundary of the brain. The interior space s is defined by the interval r ∈ [0, a), and it is characterized by the conductivity σ. The exterior to S s , nonconductive space c s is defined by r ∈ (α, ∞). The connection between the spherical and the Cartesian systems is given by

Connection Between Algorithms
In this section, we express the exterior magnetic field in spherical and ellipsoidal harmonics. The basic notation for the spectral decomposition of the Laplace operator in ellipsoidal coordinates can be found in [25], where the ellipsoidal harmonics I E m n (ρ, μ, ν) and I F m n (ρ, μ, ν) that are used in this work, as well as useful relations connecting them are explored.
The radial component of the exterior magnetic field satisfies the following relations: In spherical harmonics: In ellipsoidal harmonics: where (a 1 · a 2 · a 3 ) 1/3 a and the constants D m n , B m n denote the values of the integrals Equation (31) expressed in spherical coordinates, at a fixed point A (γ , θ , φ), γ > max{α,a 1 } gives in Cartesian form: And finally in ellipto-spherical 1 form: for each n 0,1,2 and m − 2,0,1,2.
Similarly from (32), we arrive at the relative expression in ellipsoidal coordinates at a fixed point A (γ , θ , φ), γ > max{α,a 1 } in Cartesian form: and in ellipto-spherical form: for each n 0,1,2 and m 1, 2, …, 2n + 1 Next, we identify the spherical with the ellipsoidal expansions at every point A, outside the brain. The point A, from (16), (17) and (19), is uniquely expressed in any coordinate system we use: A(r,θ ,ϕ) A(ρ,μ,ν) A(ρ,θ ε , ϕ ε ), Under this assumption and performing the necessary calculations in order to express every term appeared, in Cartesian coordinates, we arrive at the following connection relations:

Inversion MEG Problem
Let Q(r 0 ) be the moment of a dipole at the point r 0 with 0 < r 0 < c 1 . Then, the corresponding magnetic field B at the exterior domain ε satisfies [28]: Proof Recall that [36]: Using the identity Finally using the relation,

Inversion MEG Algorithm in Spherical Geometry
We assume that the magnetic field generated by a single dipole (r 0 , Q). Then the magnetic field in the exterior of the brain is provided by [1,18,26] 1 From the first 8 terms, we have Equation (59) shows that for a single dipole, since B involves r 0 × Q, the component of Q in the r direction does not contribute to the magnetic field. Thus, measurements of B yield information about the two components Q θ 0 and Q ϕ 0 . Hence, we 'lose information,' i.e., we go from two functions to a single function, as a result of Gauss's theorem.
The unique solution for the inverse EEG problem with a current dipole inside a sphere is given in (69)-(71), (79) and (80).
At this point, we use the inversion algorithm for the sphere not with its intrinsic spherical data D m n , but with the ellipsoidal data B m n as they are expressed in terms of D m n . We obtain r 0 tan ϕ 0 tan θ 0 (81)

Numerical Implementation
In the presented section, we describe what the algorithms (81)-(85), imply when we consider different values for semi-axes in order to understand the sensitivity of our analytic algorithms.
We consider an ellipsoid with semi-axes α 1 , α 2 , α 3 and a sphere with radius r α with volumes equal to that of the ellipsoid i.e V e 4π 3 α 1 α 2 α 3 V s . The conductivity is taken to be σ 0.3 ( m) −1 .
Case I We assume that there is a dipole (6) inside the ellipsoid (14) with semi-axes, α 1 9 cm a 2 6.5 cm, a 3 6 cm Therefore, the dipole source should be inside the sphere where radius is 7.054 cm but in certain directions in order to be inside the ellipsoid as well in any case.

Conclusions
For the identification of an active dipole within the brain from MEG recordings, inversion algorithms are already well known both for the spherical and the ellipsoidal model of the brain-head system. In medical practice, the developed algorithms are based on the assumption that the brain is a sphere [32]. However, the realistic model of the shape of the brain is more similar to an ellipsoid. Consequently, that generates an error on the medical practice due to the interpretation of the data, which are obtained on an ellipsoidal figure and are processed as they were obtained on a sphere. The estimation of this error is the goal of the present work. We developed the appropriate analytic formulae that connect the spherical and ellipsoidal eigenfunction expansions. Numerical calculations of the obtained results reveal that error is not significant for small eccentricities. This characteristic is clearly depicted in Figs. 1 and 2 (Tables 1,  2).
Funding Open access funding provided by HEAL-Link Greece.
Data Availability Data available within the article or its supplementary materials.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.