Helmholtz decomposition of the neuronal current for the ellipsoidal head model

In earlier work, the neuronal primary current was expressed via the Helmholtz decomposition in terms of its irrotational part characterised by a scalar function and its solenoidal part characterised by a vectorial function. Furthermore, it was shown that EEG data is affected only by the irrotational part of the current, whereas MEG data is affected by two scalar functions, namely the irrotational component and the radial part of the solenoidal vectorial function. Here, we focus on the numerical implementation of this approach on the three-layer ellipsoidal model. The parametrization of the unknown functions in terms of ellipsoidal harmonics implicitly regularizes the highly ill-posed associated inverse problems. However, despite the above parametrization of these two unknown functions in terms of ellipsoidal harmonics, the inversion matrices are highly ill-conditioned for both EEG and MEG. In order to bypass this problem, we propose an alternative approach to the inversion problem. This involves revisiting the general inversion formulas presented earlier by one of the authors and expressing them as surface integrals. By choosing a suitable parametrization for the relevant unknown functions, these surface integrals can be evaluated using a method for numerical quadrature over smooth, closed surfaces. The method uses local radial basis function interpolation for generating quadrature weights for any given node set. This gives rise to a stable linear system of equations suitable for inversion and reconstruction purposes. We illustrate the effectiveness of our approach by presenting simple reconstructions for both EEG and MEG in a setting where data are contaminated with Gaussian white noise of signal to noise ratio (SNR) of 20 dB.

Keywords: inverse problems, magnetoencephalography, EEG, MEG, numerical quadrature over surface, electroencephalography (Some figures may appear in colour only in the online journal)

Introduction
The medical significance of Electro-Magneto-encephalography, EEG-MEG, is well established, see for examples [1][2][3][4][5][6][7][8]. In particular, the high bandwidth of the EEG signal in comparison to other functional imaging modalities such as PET, SPECT and f MRI, yields a temporal resolution that offers unique insight into the study of spontaneous as well as evoked neural activity in the brain. The primary disadvantage for determining brain activity with EEG or MEG is the highly ill-posed nature of the associated spatial inverse problems: different electric currents can yield identical electric potentials measured in EEG, as well as identical magnetic fluxes measured in MEG. Actually, although it was known to Helmholtz since 1853 [9] that a current within a conductor can not be uniquely identified from the knowledge of the magnetic flux it generates in the exterior of the conductor, the first rigorous determination of the part of the current that can be obtained via MEG for the simple case of the homogeneous spherical conductor appeared in 1996 [10]. Attempts to generalize this result were made by several authors, and a variety of partial results were obtained regarding the following basic question: which part of the current affects the electric potential on the scalp and which part of the current affects the magnetic flux outside the head? In these studies, the current was assumed to be either a finite collection of dipoles or a continuous distribution of dipoles. The answer to the above question was finally obtained in [11]: let the bounded domain Ω c represent the cerebrum, which has conductivity σ c . A shell Ω f with conductivity σ f , representing the cerebrospinal fluid which surrounds the domain Ω c . The cerebrospinal fluid which is surrounded by the skull is characterized by the domain Ω b with conductivity σ b . Finally, the skull is surrounded by the scalp, which is modelled as a shell Ω s with conductivity σ s . The domain exterior to the head is denoted by Ω e , and it is assumed that Ω e is not conductive. The permeability of all domains is equal to the permeability µ of empty space. Let J p (τ ), τ ∈ Ω c , denote the primary current which is assumed to be supported within the cerebrum Ω c . This is consistent with the fact that the main source of the primary current consists of the transmembrane currents in the apical dendrites of the pyramidal cells in the cerebral cortex [5,6]. For the above situation of arbitrary geometry and arbitrary current, it is shown in [11] that the irrotational part of the current, which is characterized by the scalar function Ψ(τ ), contributes to both the electric potential on the scalp and to the magnetic field in Ω e . On the other hand, the solenoidal part of the current, which is characterized by the vectorial function A(τ ), does not affect the electric potential and furthermore only the radial component of A(τ ) affects the magnetic flux in Ω e .
It was also shown in [11] that in the case of spherical and ellipsoidal geometries, Ψ(τ ) affects the electric potential on the scalp only through its value as well as the value of ∇Ψ on the surface S c of the cerebral cortex. Similarly, the radial component of the magnetic flux outside the head is affected by the above values involving Ψ, as well as by the values of τ · A(τ ) and its gradient of S c . These results imply that in these cases the 3-dimensional nature of current distributed in a three shell head model is 'invisible' to measurements. It has been recently shown that these statements are valid for an arbitrary geometry. This fact, together with physiological considerations [5,6], have provided the motivation to study the important case that the current has support only on S c and is normal to S c . Explicit formulae expressing EEG and MEG in terms of this current will be presented elsewhere.
Several authors have focused on the dipole source problem for the ellipsoidal head model [22,23,25]. In contrast to these works, we cocentrate on the distributed source problem by decomposing the neuronal current J p into its irrotational and solenoidal components. The basic equations expressing EEG and MEG data in terms of appropriate components of the current are reviewed in section 2. In section 2.2, we show that the general formulae of inversion for both EEG and MEG which are expressed as volume integrals, can actually, be expressed as surface integrals. This paves the way for a new robust numerical implementation of the inversion equations. A brief overview of ellipsoidal geometry and ellipsoidal harmonics is presented in section 3. The inversion equations for the three layer ellipsoidal head model are given in section 4. A numerical algorithm for computing the auxiliary function v s is presented in section 5. Extensive numerical results are shown in section 6, which include numerical tests involving important terms that feature in the inversion equations derived in [11]. Furthermore, detailed numerical results and discussions are presented with regards to the properties of the inversion matrices computed from the measurement equations of [11] for both EEG and MEG. Our results are further discussed in section 7.

Head model
The different compartments of the head model are shown in figure 1. Ω c denotes the cerebrum, which is surrounded by three shells Ω f , Ω b , Ω s , denoting the cerebrospinal fluid, the skull, and the scalp. Their conductivities are respectively denoted by σ c , σ f , σ b and σ s . The spaces Ω c , Ω f , Ω b and Ω s are bounded by the surfaces S c , (S c , S f ), (S f , S b ) and (S b , S s ), respectively. Table 1 presents the conductivity values of the head model as documented in [12][13][14].
The physiology of the head model is accurately characterized by the four layer compartments shown in figure 1. The distributed inverse source results derived in [11] consider the four layer head model. It can be observed from table 1 that Ω f (CSF), has a higher conductivity than the remaining compartments but it also has a very small thickness. The detailed analysis of [15] shows that the brain-CSF interface has a negligible effect in the forward model. For this reason, in the case of our numerical examples, we ignore Ω f (CSF) and restrict our analysis to the three layer head model involving the compartments Ω c , Ω b and Ω s .

Auxiliary functions v j (r, τ )
It was shown in [11] that for a given geometry, the functions v j (r, τ ), j = c, f , b, s are defined via the following boundary value problem: Equations (1)-(4) are independent of the current J p (τ ) and depend only on the geometry and on the conductivities σ c , σ b , σ f and σ s .
It is shown in [11] that the functions v j (r, τ ) can be related to the functions u j (r, τ ) with unit Volts, r ∈ Ω j , τ ∈ Ω c , j ∈ {c, f , b, s}. They are defined in terms of a single dipole with moment Q(τ ) with unit Coulomb-meter, located at the position vector τ via the following equations: The functions u j and v j are related by the equation

Formulae for a current with three-dimensional support
Let v s (r, τ ) be defined by equations (1)-(4) and suppose that the primary current J p (τ ), τ ∈ Ω c , is continously distributed in the cerebral cortex. It is shown in [11] that the electric potential u s (r) on the scalp is given by It was assumed in [11] that the support of the primary current is a closed set inside the open set Ω c , and therefore it was assumed that there exists an > 0 such that the shell S c− of thickness interior to Ω c is free of neuronal sources. Thus, it was assumed that J p (τ ) vanishes on S c . Under the assumption that J p (τ ) has sufficient smoothness, we can use the Helmholtz decomposition. The Helmholtz decomposition of an arbitrary vectorial function implies that J p can be represented in the form where A(τ ) satisfies the constraint ∇ · A(τ ) = 0. Here, A(τ ) and Ψ(τ ) also have sufficient smoothness. The constraint (∇ · A(τ ) = 0) implies that J p (τ ) involves three arbitrary scalar functions, namely the scalar function Ψ(τ ) and the two independent scalar functions characterising A(τ ). Replacing in the rhs of (10) the divergence of J p (τ ) with the Laplacian of Ψ(τ ), we find It is also shown in [11] that the radial component of the magnetic flux in the space exterior to the head, which is denoted by Ω e , is given by

Ωc
(r · H(r, τ ))(∆ τ Ψ(τ ))dV(τ ), r ∈ Ω e (13) where H(τ ) is explicitly defined in [11] in terms of v j , j = c, f , b, s. Green's identity states that if u and v are sufficiently smooth, then Using this identity with we can map the volume integral appearing in the rhs of (12) to a surface integral. Indeed, the equations defining the functions v j (r, τ ), j = c, f , b, s, remain invariant under the interchange of r with respect to τ . Thus, v s is also harmonic with respect to τ . Thus, using equation (14) with the identifications (15), equation (12) becomes Under the assumption that J p (τ ) vanishes on S c , which was justified earlier, equation (16) becomes (17) Thus, in this case u(r) is affected only by the values of Ψ and of ∇ τ Ψ on the surface of S c .
For the particular case of spherical geometry, r · H(r, τ ) = 0 in the second term on the right hand side of equation (13), see [18]. For arbitrary head models, the contribution of this second term is negligible. So, equation (13) reduces to Employing equation (14) with in equation (18) we arrive at the following surface inversion formula for MEG: where A τ (τ ) := τ · A(τ ). Thus, r · B(r) is affected only by the values of A τ (τ ) and of ∇ τ A τ (τ ) on the surface of S c .
In the particular cases of spherical and ellipsoidal geometries, equations (17) and (19) were actually already obtained in [11] (see equations (4.6) and (5.14) of [11]). However, it was not realized in [11] that equations (17) and (19) are actually valid for an arbitrary geometry.

Ellipsoidal geometry
Ellipsoidal shells provide a good approximation for a human brain [11,[22][23][24][25][26]. It is well known that ellipsoidal harmonics are far more complicated than spherical harmonics [16]. However, the comprehensive paper of Bardhan and Knepley [27] facilitates the computational implementation of ellipsoidal harmonics for a large class of problems ranging from molecules to solar systems. The above paper elucidates the various challenges occurring in ellipsoidal harmonics, including the sign ambiguities associated with different coordinate transformations. It also discusses limitations associated with certain numerical implementations, and it includes readily available open source codes for both the Matlab and Python platforms.
In what follows we adopt the notations of [27]. An ellipsoidal surface satisfies the equation where the constants a, b, c called the ellipsoidal semi-axes, satisfy the inequalities The term ellipsoidal is used for a tri-axial ellipsoidal (three unequal semi-axes).

Ellipsoid coordinates and separation of variables
In the ellipsoid coordinate system, a point r = (x, y, z) in Cartesian coordinates is written as (λ, µ, ν) in ellipsoidal coordinates. Each ellipsoidal coordinate is a root of the following equation which is a cubic algebraic equation in the variable s 2 : with The squares of the ellipsoidal coordinates are in the following ranges: Points on the surface of the ellipsoid with semi-axes a, b and c, satisfy the equation λ = a, where λ ∈ [k, +∞). However, the authors of [27] do not enforce the non-negativity assumption as this poses a problem for inverse coordinate transformation (this issue is discussed in sections 3 and 5 of [27]). For a given point r = (x, y, z), the magnitude of the ellipsoidal coordinates (λ, µ, ν) can be computed via the following equations: where The Cartesian coordinates can be expressed in terms of the ellipsoidal coordinates via the equations

The Lamé equation and its solutions
In ellipsoidal harmonics, the Laplace equation separates in a symmetric way, so that the solution of each coordinate satisfies the same differential equation, which is called Lamé's equation: where The functions F p n and E p n depends on h and k, thus they are invariant for a set of co-focal ellipsoidal surfaces.

Ellipsoidal harmonics
For a given degree n and order p, the interior and exterior ellipsoidal harmonics are defined by and respectively. The surface ellipsoid harmonics are defined as The surface harmonics satisfy the orthogonality condition where the normalization constants γ p n are given by The Coulomb potential at r due to a unit charge at τ with |r| > |τ |, can be expanded in terms of ellipsoidal harmonics via the equation The normal derivative at the ellipsoid surface defined by λ = a is computed via the equation If n is the unit normal to the ellipsoid surface λ = a, then Equations (34) and (36) are verified numerically in section 6.

Measurement equations for the three layer ellipsoidal head model
The objective in this section is to present equations (12) and (13) for the particular case of the three layer ellipsoidal head geometry (ignoring the CSF layer). The ellipsoidal harmonics are computed using the open source code provided by [27]. The surfaces S c , S b , S s shown in figure 1 are now co-focal surfaces with the following characteristics: The surface of the cerebrum S c is defined by the equation Using, the results of [11,23,25,27] together with extensive numerical testing, we have obtained the following numerically stable representation for the auxiliary function v s (r, τ ): where C p n are geometry dependent coefficients, r = (λ r , µ r , ν r ) is the position vector of the electrode and γ p n is given by equation (33). In other words, in order to obtain a numerically stable approach for computing the expansion coefficients of v s , it is convenient to express these coefficients in the form C p n γ p n . A numerically robust approach for estimating the coefficients C m l using a BEM solver will be discussed in section 5.
It is shown in [11] that the function H(r, τ ) defined in equation (2.17) of [11] can be expressed in the following form: where the Cartesian components of the constant vector H m,l n depend on the conductivities and on the parameters of the co-focal surfaces 3}}. An explicit formula for H m,l n is given in the appendix of [11], however, H m,l n can be computed directly, in analogy with the approach for computing C m l . Fortunately, the contribution of the second term of (13) is much smaller than the contribution of the first term of the rhs of (13), thus in practice the second term can be neglected.
For τ ∈ Ω c , the function Ψ(τ ) can be expanded in the form It is shown in [11] that equation (12) yields the following equation for the case of the ellipsoidal geometry: Similarly, equation (13) yields the equation

Estimating the geometry dependent coefficients C m l
The coefficients C m l appearing in (38) have been derived in [23]. However, unlike the analogous expression in the case of spherical geometry derived in [19], the expressions derived in [23] are far more complex. We propose to estimate the relevant geometry-dependent coefficients C m l , using data generated from a BEM solver: in the case of an ellipsoidal head model, using equations (9) and (38) the electric potential u s (r, τ , Q(τ ), C m l ) due to a dipole source Q(τ ) is given by where ∇ τ E m l (τ ) can be estimated accurately, using a finite difference approximation (numerical results regarding the accuracy of a finite difference approximation of ∇ τ E m l (τ ) are given in section 6, see figure 4). There exist several freely available numerical solvers for the boundary value problem described by equations (5)- (8) for an arbitrary head model, see for example [20,21]. The steps for a numerical estimate of C m l are as follows: (i) Choose N source-observation pairs (τ , r) of randomly oriented dipoles Q(τ ).
(ii) Employ a BEM/FEM solver to compute the potentials for each source observation pair (τ , r). Here, we use the freely available openMEEG solver openMEEG [20] to obtain these potentials and denote them by {ũ i : i = 1, ...N}. These solutions feature as data in a minimization algorithm. (iii) It is clear from equation (45), that the relationship between u s (r, τ , Q, C m l ) and C m l is linear. The least squares estimate of C m l is obtained by the minimizing function where u s (r, τ , Q, C m l ) denotes the parametric form of equation (45) and {ũ i : i = 1, ..., N} are the data generated in step (ii). The estimated geometry-dependent coefficients using the proposed numerical procedure are denoted by C m l to differentiate them from the exact coefficients denoted by C m l . (iv) The above procedure can be cross validated as follows: Choose a new set of sourceobservation pairs (τ , r) of randomly oriented dipoles Q(τ ). Employ a BEM/FEM solver to compute the potential for each source observation pair (τ , r) to generate the data set {u i : 1 i T}. Then, employing the estimated C m l (step 3) in equation (45) to generate estimated data u s (r, τ , Q,C m l ) and compare them with the cross validation data. Results of such comparison are shown in section 6, see figure 5.

Numerical results
Numerical tests relevant to the implementation of ellipsoidal harmonics employed in this study are presented in [27]. Here we present additional numerical results for several important terms which feature in the inversion formulae of [11] with the relevant root mean square errors (RMSE). Following the recommendations of the authors of [27] regarding numerical stability, we have restricted the order of the ellipsoidal harmonics to n 12.
The parameters of the three layer ellipsoidal head model used in this study are given in table 3. These parameters are obtained by fitting ellipsoidal shells to a realistic head model.

Numerical verification of (34) and (36)
In this section we compare numerically the left and right hand sides of equations (34) and (36). This provides comparisons between analytic expressions and their corresponding ellipsoidal harmonic expansions, where the maximum degree is n max = 12. These comparisons provide important insight into the convergence properties of ellipsoidal harmonics and the accuracy of the implementation described in [27]. We consider two different cases, relevant to EEG and MEG.
In the first case, we place the observation point on the surface of the scalp, i.e. r ∈ S s . This is relevant to EEG, where the sensor array is positioned on the scalp. Results are shown in  The second test case, is identical to the first case except that the observation point is placed outside the head, i.e. r ∈ Ω e . This is relevant to MEG, where the sensor array is outside the head. Results are shown in A comparison of the absolute errors of figures 2 and 3 (subplots (c) and (d)) suggest the numerical stability in the implementation of ellipsoidal harmonics as outlined in [27] is sensitive to the distance between the source point τ ∈ S c and the observation point r, i.e. d = |r − τ |. More precisely, the errors between the left and right hand sides of equations (34) and (36) increase substantially as the distance between the source-observation pair (τ , r) decreases. Moreover, in tests involving increasing the number of terms in the ellipsoidal harmonic expansion for the right hand sides of equations (34) and (36), we have found that the numerics become very unstable for n max > 12, consistent with discussions in [27].

The computation of C m l
Equation (45) requires the computation of the gradient ∇ τ E p n (τ ), which we choose to estimate using the simple finite difference technique. A reasonable test in assessing this estimate (∇ τ E p n (τ )) is to consider again a simple test involving equation (34): we consider the scalar product of the gradient of equation (34) with a dipole Q(τ ). We place 462 dipoles Q(τ ) in the cerebrum, τ ∈ Ω c , and choose an arbitrary observation vector on the scalp surface r ∈ S s . The test involves comparing Q · ∇ τ 1 |r−τ | (which is straightforward to compute) with its corresponding ellipsoidal harmonics approximation given below The finite difference technique is used to estimate ∇ τ E p n (τ ) occurring in the right hand side of (47). The results are shown in figure 4.

Cross validation of C m l
In section 5, a numerical procedure was outlined for computing the geometry-dependent coefficients C m l appearing in equation (38). These coefficients are central to the forward problem. The estimated coefficients C m l are then used to generate a new data set which is compared with data generated from a BEM solver [20]. This approach is referred to as cross validation. A total of 9700 source-observation pairs were used during the estimation step, i.e. N = 9700 in equation (46). The geometric and mesh parameters used in the BEM computations are given in table 3.
The results of this analysis are shown in figure 5. Subplot (a) depicts the estimated coefficients log[C m l ] presented in sequential order. For example, on the x-axis, index 10, refers to the coefficient C 1 4 . Subplot (b) depicts the comparison between the estimated data using equation (45) and data generated using a BEM solver, on a new set of source-observation pairs. Subplot (c) shows the absolute error, i.e. Error(i) := |u i − u s (r i , τ i , Q,C m l )|, where C m l have been used to generate the estimated data.

EEG inversion matrix
In this section, we present reconstructions for EEG using synthetic data. Recall the EEG inversion formula given by equation (43). For reasons related to the sign ambiguity issue and discussed in [27], we rewrite the parametrization of equation (41) as where the functions ψ p n (λ) have a λ dependence. For the particular case of the spherical head model, for both EEG and MEG a procedure was outlined in [17] to establish the radial dependence via a minimization procedure. The expansion of data on an ellipsoidal surface fixes the angular dependence, thus only the radial dependence is not determined. Here, for simplicity we assume that ψ p n (λ) takes the form ψ p n (λ) = λa p n . In future studies this dependence will be fixed via the minimization of L 2 norm. Equation (48) is employed in the numerical computations due to problems associated with sign ambiguity in the implementation of ellipsoidal harmonics discussed in [27]. On the surface S c , λ = a 3 , and {a p n } are the unknown coefficients linearly related to data. Equation (43) gives rise to the following system of linear equations: where u is the data (from EEG electrodes), E ∈ R Ns×Np is the ellipsoidal inversion matrix constructed exactly as stated by equation (43), and α = {a p n } ∈ R Np is the vector of unknown coefficients to be estimated from the data. Most electrode caps have between 60 to 128  electrodes. However, for the sake of a more rigorous analysis and understanding of the properties of the ellipsoidal inversion matrix, we have assumed a sensor array of 312 electrodes distributed evenly around the scalp, i.e. N s = 312. Indeed, this is a theoretical study aiming to improve our understanding in preparation for the practical implementations. The position vectors of the electrodes on the surface S s are simply a subset of the nodes of the surface S s . Surface S s has N = 3121 nodes and one electrode has been placed for every tenth node of this mesh. This is an ideal scenario. In practice, it is physically impossible to place one electrode at every tenth node. Despite this ideal situation, the results will show that this problem is still highly ill-conditioned. Moreover, we have restricted the dimension of the vector of unknown coefficients a p n to be N p = 100, i.e. α ∈ R 100 . So, we are intentionally considering the case that we have more data than unknowns.
Our first numerical test involves a very simple test of generating a vector of random coefficients α = {a p n } to obtain a data set via equation (49), u = Eα. We then obtain an   For example on the x-axis, index 5 refers to (l = 2, m = 0) and index 25 on the x-axis, refers to (l = 4, m = 8), i.e. C 8 4 . In these computations, n max = 12. Subplot (b) depicts the comparison between test data and estimated data u s (r, τ , Q,C m l ). Subplot(c) depicts the absolute error between test and estimated data depicted in Subplot (b) with a RMSE = 1.7285. n and order p. So, the contribution to the data from higher degree ellipsoidal harmonic terms is very small. For example, on average, the matrix entry in E * ,j multiplying the coefficient a 1 1 , is 3 × 10 12 larger than the entry multiplying a 8 4 , i.e. E k,3 E k,25 ≈ 3 × 10 12 , ∀k N s . This fact poses a severe restriction to the class of functions that can be represented by the expansion of equation (41). Figure 6(c) depicts the logarithmic plot of the singular values of the ellipsoidal inversion matrix E which is consistent with the conclusion drawn from figure 6(b). It shows that the singular values decay quickly. The inversion matrix E has a rank of 23 with a condition number matrix κ(E) = 4.1 × 10 62 . Despite the fact that the above problem is overdetermined, the matrix E is severely ill-conditioned and should not be employed directly for inversion or reconstruction purposes. In the case of the spherical model problem, the inverse problem was numerically solved via the technique of reproducing kernels [18]. However, deriving such kernels for the ellipsoid head model is challenging and is work in progress. It is worth noting that equation (43) may still be employed to generate data for a given restricted class of Ψ(τ ) functions.
In summary: first, the expansion in terms of ellipsoidal harmonics is suitable only for a very restricted class of functions, and second, this expansion yields an inversion matrix with a very high condition number.

EEG reconstructions
Due to the above difficulty we will not solve the inverse problem using the ellipsoidal harmonics and equation (49). Instead, we construct an inversion matrix using the surface inversion equation (17), and a numerical quadrature technique involving local radial basis functions as outlined in [29]. This important paper outlines a numerical scheme for the integration of functions on any smooth, closed, triangulated surface. This technique only requires function values of the integrands at the nodes of the triangulated mesh. Moreover, the paper is complemented with an open source Matlab implementation which includes detailed numerical tests. In this setting, we choose to parametrize Ψ(τ ) using multi-quadratic radial basis functions: where {β j } are the radial basis coefficients and c is the shape parameter. The shape parameter should be estimated from the data. However, to avoid high computational costs, we have abitrarily set this parameter to be c = 10 −2 . In order to avoid the inverse crime [28], we generate the data using the ellipsoidal harmonic inversion equation (49). This restricts the class of functions analysed here, but the inversion method introduced in this section is valid for any smooth function. The construction of the inversion matrix involves employing equations (50) in (17). The surface integral is evaluated using the quadrature integration method outlined in [29]. We consider the case that data is contaminated with additive white Gaussian noise (WGN) with a signal to noise ratio of 20 dB (SNR = 20 dB). This leads to the following system of linear equations where L ∈ R Ns×N . Here N s = 312 is the number of electrodes and N = 3124 and is the number of nodes in the surface mesh of the cerebrum (see table 3 for details of the surface meshes of the head ellipsoidal head model).
Recall that figure 6(b) depicts log[|E ij |], which decays quite rapidly with increasing degree n, and order p. Moreover, figures 6(a) and (c) suggest that only 25 coefficients can be recovered. So, in this simple test case we consider that data is generated from only the coefficient a 8 4 , which corresponds to column entry 25 that multiplies the matrix entries E k, 25 , for all sensor indices k. All other coefficients are set to zero, i.e. {a p n = 0 : n = 4, p = 8}. The reconstruction steps are outlined below:   Figure 11. This figure depicts a reconstruction test case on an ellipsoidal head model for MEG. The signal to noise ratio is SNR = 20 dB. Subplot (a) depicts the comparison of A τ (τ ) and its corresponding reconstruction Â τ (τ ). Blue solid line with circle is A τ (τ ) and the red solid line with cross is Â τ (τ ). Subplot (b) depicts the absolute error E( j) := |A τ (τ j ) −Â τ (τ j )| for every node index with a RMSE = 0.2206.

MEG inversion matrix
In this section, we present reconstructions using synthetic data for MEG. Recall, from equation (42) that τ A τ (τ ) can be parametrized in the form below We follow an identical procedure to the one outlined in section 6.4 to study the properties of the MEG inversion matrix. A typical MEG machine has generally 102 magnetometers. Here we assume an array of 518 magnetometers distributed evenly outside the head. On average, we placed one magnetometer r ∈ Ω e for every six nodes of the surface S s . Despite employing far more magnetometers than a MEG machine, the results will show that the relevant inversion matrix is highly ill-conditioned. The inversion formula used for this reconstruction is given by equation (44). The contribution from the second term in equation (44) involving r · H m,l n is neglegible and will be omitted in this analysis, see [18]. For simplicity, we assume that the functions f p n (λ) take the form f p n (λ) = λa p n . On the surface S c , λ = a 3 and {a p n } are the unknown coefficients linearly related to data, i.e. to the left hand side of equation (44). Equation (44) gives rise to the following system of linear equations: where b is the data (from MEG electrodes), M ∈ R Ns×Np is the MEG inversion matrix constructed from equation (44), neglecting the contribution from the second term. The vector η = {a p n } ∈ R Np is the vector of unknown coefficients to be estimated from data. Most MEG sensor arrays are between 108 to 306 channels. We have restricted the dimension of the vector of unknown coefficients η = {a p n } to be N p = 100, i.e. β ∈ R 100 . So, we are considering the case that we have more data than unknowns.
Similar to the case of EEG of section 6.4, we generate a vector of random coefficients η = {a p n } to obtain a data set via equation (53), b = Mη. We proceed to obtain the estimate η via the same equation, η = M −1 b. In this test, we have not added noise to the data. The results are shown in figure 9. Figure 9(a) shows the comparison of the randomly generated vector β and its corresponding estimate η. It can be observed that we can only recover at most around 22 coefficients out of 100. This corresponds to the degree n and order p of the ellipsoidal harmonics E p n being in the set {(n, p) : 0 < n 4, 0 p 2n}. Moreover, {(n = 0, p = 0), (n = 1, p = 0)} are also excluded from this set. This result can be explained by figure 9(b), which shows the logarithmic plot of all entries of the matrix, i.e. {log[|M i,j |] : i N s , j N}. This figure shows that for all sensors (rows of M), the entries multiplying the relevant coefficients {a p n } decay very quickly. So, the contribution of higher degree ellipsoidal harmonic terms are very small on the data. For example, on average, the matrix entry in E * ,j multiplying the coefficient a 1 1 is 2.8 × 10 8 larger than the entry multiplying a 5 4 , i.e. M k,3 M k,22 ≈ 2.8 × 10 8 , ∀k N s . This limits the class of functions that can be represented by the expansion of equation (42). Figure 9(c) shows the logarithmic plot of the singular values of the matrix M which is consistent with the conclusion drawn from figure 9(b): it shows that the singular values decay quickly and the rank of this matrix is 22 with a condition number of κ(M) = 2.38 × 10 29 . As in the case of EEG, this matrix is highly ill-conditioned and difficult to invert.
The reconstruction steps are the same steps taken for the case of EEG as discussed in section 6.5. However, here we set a 3 4 = 1 and the remaining a p n = 0. The data is generated via equation (53). For the inversion step, an inversion matrix is constructed using the surface inversion equation (19) and a numerical quadrature technique involving local radial basis functions as outlined in [29]. The reconstructions for MEG on an ellipsoidal head model are depicted in figure 10. Figure 11(a) depicts the comparison between the test function and the reconstructed function for every node on the surface of the cortex and figure 11(b) depicts the corresponding absolute error.

Conclusion
The modeling of dipole sources for an ellipsoidal head geometry has been covered rigorously in the important studies of [22,23,25]. Here, we analyse the case that J p (τ ) is continously distributed in Ω c [11]. In the particular case of the ellipsoidal head model, the assumption that ellipsoidal harmonics represent a suitable basis for the expansion of the unknown functions Ψ(τ ), A τ (τ ) may be interpreted as implicit regularization and an attempt to mitigate the ill-posedness of the associated inverse problems. However, we show that the inversion matrices constructed from equations (43) and (44) are highly ill-conditioned and unsuitable for reconstructions. This is the case even for a highly overdetermined system of equations (more data than unknowns). Robust inversion techniques such as the reproducing kernel method were successful for the EEG and MEG inversion problem in the case of the spherical head model [18]. However, deriving these kernels for the ellipsoidal geometry is more challenging and is currently under study. The above difficulties can be overcome by using an alternative approach to the inversion equation. This involves reformulating the volume integral equations of (12) and (13) as the surface integrals of (17) and (19). In this setting, Ψ(τ ) and A τ (τ ) were parametrized in terms of multiquadratic radial basis functions. The surface integrals were computed using the powerful numerical quadrature technique outlined in [29]. This gave rise to inversion matrices suitable for reconstructions. In order to avoid the 'inverse crime', data were generated using the ellipsoidal harmonic expansions and equations (43) and (44). Reconstructions shown in figures 7 and 10 suggest that the above approach based on the inversion equations (17) and (19) is a new, robust, reconstruction approach to the distributed inverse source problem in EEG and MEG.