Effects of head modeling errors on the spatial frequency representation of MEG

It is common practice in MEG and EEG modeling to express the total current as a superposition of two components, i.e.

$$\begin{eqnarray}&&{\bf{J}}\left({\bf{r}}^{\prime} \right)={{\bf{J}}}^{P}\left({\bf{r}}^{\prime} \right)+{{\bf{J}}}^{v}\left({\bf{r}}^{\prime} \right),\end{eqnarray} \tag{ 1 }$$

where ${{\bf{J}}}^{P}\left({\bf{r}}^{\prime} \right)$ represents the physiologically interesting primary current and ${{\bf{J}}}^{v}\left({\bf{r}}^{\prime} \right)$ is the associated passive volume current component that completes the loop of electric current in the brain tissue. On the basis of the Coulomb force, the passive volume current is of the form

$$\begin{eqnarray}&&{{\bf{J}}}^{v}\left({\bf{r}}^{\prime} \right)=\sigma {\bf{E}}=-\sigma {\rm{\nabla }}V,\end{eqnarray} \tag{ 2 }$$

where σ, E, and V are the electric conductivity, the electric field, and the electric potential, respectively. The second equality in equation (2) is valid when the time derivative of the magnetic field B is insignificant, which is the case in MEG and EEG (Hämäläinen et al 1993). If we assume a piecewise homogeneous conductor head model with N_S conductivity boundary surfaces, the magnetic field due to ${\bf{J}}\left({\bf{r}}^{\prime} \right)$ is given by Geselowitz' formula (Geselowitz 1967, Hämäläinen et al 1993)

$$\begin{eqnarray}&&{\bf{B}}\left({\bf{r}}\right)={{\bf{B}}}_{0}\left({\bf{r}}\right)+\displaystyle \frac{{\mu }_{0}}{4\pi }\displaystyle \sum _{l=1}^{{N}_{S}}\left({\sigma }_{l}^{-}-{\sigma }_{l}^{+}\right){\int }_{{S}_{l}^{\prime} }V\left({\bf{r}}^{\prime} \right)\displaystyle \frac{{\bf{r}}-{\bf{r}}^{\prime} }{{\left|{\bf{r}}-{\bf{r}}^{\prime} \right|}^{3}}\times d{\bf{S}}{{\prime} }_{l},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{{\bf{B}}}_{0}\left({\bf{r}}\right)=\displaystyle \frac{{\mu }_{0}}{4\pi }{\int }_{v^{\prime} }{{\bf{J}}}^{P}\left({\bf{r}}^{\prime} \right)\times \displaystyle \frac{{\bf{r}}-{\bf{r}}^{\prime} }{{\left|{\bf{r}}-{\bf{r}}^{\prime} \right|}^{3}}\ {dv}^{\prime} \end{eqnarray} \tag{ 4 }$$

is the contribution by primary currents, $v^{\prime}$ is the total head volume, and ${\sigma }_{l}^{-}$ and ${\sigma }_{l}^{+}$ denote conductivities of the inner and outer regions relative to S_l . The second term on the right hand side of (3) is the contribution by volume currents, which we will denote as B_vol . We need the electric potential V in order to calculate B_vol and hence the total magnetic field. Below we review the expressions to obtain discrete values of V on the boundaries that we may thus use to approximate and discretize B.

Here, we give an overview of the CC BEM approach. The LC requires only a slight modification, as will be pointed out. More comprehensive reviews may be found in (de Munck 1992, Ferguson et al 1994, Phillips et al 2006, Stenroos et al 2007, Stenroos 2009).

The electric potential V(r) for a field point r on the k^th surface S_k is given intrinsically by Hämäläinen et al (1993), Phillips et al (2006)

$$\begin{eqnarray}&&V\left({\bf{r}}\right)={V}_{\infty }\left({\bf{r}}\right)-\displaystyle \frac{1}{2\pi }\displaystyle \sum _{l=1}^{{N}_{S}}\displaystyle \frac{{\sigma }_{l}^{-}-{\sigma }_{l}^{+}}{{\sigma }_{k}^{-}+{\sigma }_{k}^{+}}{\int }_{{S}_{l}^{\prime} }V\left({\bf{r}}^{\prime} \right)\displaystyle \frac{{\bf{r}}-{\bf{r}}^{\prime} }{{\left|{\bf{r}}-{\bf{r}}^{\prime} \right|}^{3}}\cdot d{\bf{S}}{{\prime} }_{l},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&{V}_{\infty }\left({\bf{r}}\right)=\displaystyle \frac{1}{2\pi \left({\sigma }_{k}^{-}+{\sigma }_{k}^{+}\right)}{\int }_{v^{\prime} }{{\bf{J}}}^{P}\left({\bf{r}}^{\prime} \right)\cdot \displaystyle \frac{{\bf{r}}-{\bf{r}}^{\prime} }{{\left|{\bf{r}}-{\bf{r}}^{\prime} \right|}^{3}}\ {dv}^{\prime} .\end{eqnarray} \tag{ 6 }$$

Like in the B field case, equation (6) describes the primary current contribution to the potential whereas the surface integral is the volume current contribution. The primary current contribution is easily obtained if we have prescribed J^P . In particular, if we let J^P be a collection of current dipoles with positions represented by delta functions, the integral collapses to a simple form. For N dipoles ${{\bf{J}}}_{n}^{P}={{\bf{Q}}}_{n}\delta ({\bf{r}}-{{\bf{r}}}_{n})$ where n = 1, ..., N, we have

$$\begin{eqnarray}&&{V}_{\infty }\left({\bf{r}}\right)=\displaystyle \frac{1}{2\pi \left({\sigma }_{k}^{-}+{\sigma }_{k}^{+}\right)}\displaystyle \sum _{n=1}^{N}{{\bf{Q}}}_{n}\cdot \displaystyle \frac{{\bf{r}}-{{\bf{r}}}_{n}}{{\left|{\bf{r}}-{{\bf{r}}}_{n}\right|}^{3}}.\end{eqnarray} \tag{ 7 }$$

As for the volume current contribution, we may discretize each surface S_l into N_l triangles. Then, (5) can be written as

$$\begin{eqnarray}&&V\left({\bf{r}}\right)={V}_{\infty }\left({\bf{r}}\right)-\displaystyle \frac{1}{2\pi }\displaystyle \sum _{l=1}^{{N}_{S}}\displaystyle \frac{{\sigma }_{l}^{-}-{\sigma }_{l}^{+}}{{\sigma }_{k}^{-}+{\sigma }_{k}^{+}}\displaystyle \sum _{m=1}^{{N}_{l}}{\int }_{{{\rm{\Delta }}}_{l}^{m}}V\left({\bf{r}}^{\prime} \right)\displaystyle \frac{{\bf{r}}-{\bf{r}}^{\prime} }{{\left|{\bf{r}}-{\bf{r}}^{\prime} \right|}^{3}}\cdot d{{\bf{S}}}_{{{\rm{\Delta }}}_{l}^{m}}^{\prime} ,\end{eqnarray} \tag{ 8 }$$

where ${{\rm{\Delta }}}_{l}^{m}$ is the mth triangle of surface S_l .

If we have chosen a large enough N_l such that the triangle areas are small, we may reasonably make some assumptions about the behavior of the potential V within each triangle. In turn, V can be estimated with some basis and weight functions defined relative to parameters of relevant triangles. The potential V can then be explicitly defined and solved for. Many such approximations exist, including using constant or linear basis with collocation or Galerkin weighting (Stenroos 2009). Higher-degree basis functions have been considered as well (Gençer and Tanzer 1999).

In this paper, we use BEM as a tool to calculate the errors of volume current contributions due to closer sensors. Since we are not actually concerned with the accuracy between different approximation methods but rather the overall behaviour of the resulting signal due to head geometry errors, the most straightforward approximations which are the CC and LC approaches suffice for our purposes. We present the CC case in the next section which assumes a constant potential within each triangle. The LC approach approximates the potential within each triangle as a linear function via an interpolation from the potentials at the three vertices; see de Munck (1992), Stenroos et al (2007), Stenroos (2009). The simplicity of the CC approach allows us to outline an analytical method one may use to find errors for small boundary perturbations as well (see the supplementary material).

First, we may assume that triangles are small enough such that the potential is constant in each triangle, i.e. $V({\bf{r}}^{\prime} )\approx V({{\bf{c}}}_{l}^{m})$ when ${\bf{r}}^{\prime} \in {{\rm{\Delta }}}_{l}^{m}$, where ${{\bf{c}}}_{l}^{m}$ is the centroid of ${{\rm{\Delta }}}_{l}^{m}$. This allows us to pull the potential term out of the integral, and we get

$$\begin{eqnarray}&&{\int }_{{{\rm{\Delta }}}_{l}^{m}}V\left({\bf{r}}^{\prime} \right)\displaystyle \frac{{\bf{r}}-{\bf{r}}^{\prime} }{{\left|{\bf{r}}-{\bf{r}}^{\prime} \right|}^{3}}\cdot d{{\bf{S}}}_{{{\rm{\Delta }}}_{l}^{m}}\approx V\left({{\bf{c}}}_{l}^{m}\right){\int }_{{{\rm{\Delta }}}_{l}^{m}}\displaystyle \frac{{\bf{r}}-{\bf{r}}^{\prime} }{{\left|{\bf{r}}-{\bf{r}}^{\prime} \right|}^{3}}\cdot d{{\bf{S}}}_{{{\rm{\Delta }}}_{l}^{m}}^{\prime} .\end{eqnarray} \tag{ 9 }$$

Notice that the integral on the right hand side is now simply the solid angle spanned by the triangle ${{\rm{\Delta }}}_{l}^{m}$ from the observation point r; let us denote it as ${{\rm{\Omega }}}_{l}^{m}\left({\bf{r}}\right)$. If we let r_i0 ≡ r_i − r, i = 1, 2, 3 be the three vertices of the triangle relative to r, and let r_i0 be their lengths, we may equivalently express each ${{\rm{\Omega }}}_{l}^{m}$ as (Oosterom and Strakee 1983)

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{l}^{m}=2\arctan \left[\displaystyle \frac{{{\bf{r}}}_{10}\cdot \left({{\bf{r}}}_{20}\times {{\bf{r}}}_{30}\right)}{{r}_{10}{r}_{20}{r}_{30}+\left({{\bf{r}}}_{10}\cdot {{\bf{r}}}_{20}\right){r}_{30}+\left({{\bf{r}}}_{30}\cdot {{\bf{r}}}_{10}\right){r}_{20}+\left({{\bf{r}}}_{20}\cdot {{\bf{r}}}_{30}\right){r}_{10}}\right].\end{eqnarray} \tag{ 10 }$$

Let r coincide with centroids of the triangles as well. Then, all the potential terms of (5) are discretized at the same locations and it can now be compactly written as

$$\begin{eqnarray}&&V\left({{\bf{c}}}_{k}^{i}\right)={V}_{\infty }\left({{\bf{c}}}_{k}^{i}\right)-\displaystyle \frac{1}{2\pi }\displaystyle \sum _{l=1}^{{N}_{S}}\displaystyle \sum _{m=1}^{{N}_{l}}\displaystyle \frac{{\sigma }_{l}^{-}-{\sigma }_{l}^{+}}{{\sigma }_{k}^{-}+{\sigma }_{k}^{+}}V\left({{\bf{c}}}_{l}^{m}\right){{\rm{\Omega }}}_{l}^{m}\left({{\bf{c}}}_{k}^{i}\right).\end{eqnarray} \tag{ 11 }$$

In matrix form, this looks like

$$\begin{eqnarray}&&\left[\begin{array}{c}{{\bf{V}}}_{1}\\ \vdots \\ {{\bf{V}}}_{{N}_{S}}\end{array}\right]=\left[\begin{array}{c}{{\bf{V}}}_{\infty ,1}\\ \vdots \\ {{\bf{V}}}_{\infty ,{N}_{S}}\end{array}\right]+\left[\begin{array}{ccc}{{\bf{G}}}_{\mathrm{1,1}} & \cdots & {{\bf{G}}}_{1,{N}_{S}}\\ \vdots & \ddots & \vdots \\ {{\bf{G}}}_{{N}_{S},1} & \cdots & {{\bf{G}}}_{{N}_{S},{N}_{S}}\end{array}\right]\left[\begin{array}{c}{{\bf{V}}}_{1}\\ \vdots \\ {{\bf{V}}}_{{N}_{S}}\end{array}\right],\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&{G}_{k,l}^{i,m}=-\displaystyle \frac{1}{2\pi }\displaystyle \frac{{\sigma }_{l}^{-}-{\sigma }_{l}^{+}}{{\sigma }_{k}^{-}+{\sigma }_{k}^{+}}{{\rm{\Omega }}}_{l}^{m}\left({{\bf{c}}}_{k}^{i}\right).\end{eqnarray} \tag{ 13 }$$

Note that the G matrix is dependent only on the geometry and conductivities of the conductor, and it is also the only term that depends on the boundary geometry. This means that for each different source configuration in the same head model, G needs to be calculated only once, whereas the primary current contribution V_∞ needs to be recalculated.

If we write the matrix equation (12) as V = V_∞ + GV, then we have

$$\begin{eqnarray}&&\left({\mathbb{I}}-{\bf{G}}\right){\bf{V}}={{\bf{V}}}_{\infty },\end{eqnarray} \tag{ 14 }$$

where ${\mathbb{I}}$ is the identity matrix. It is seemingly straightforward to solve for V by taking the inverse of $({\mathbb{I}}-{\bf{G}})$, but $({\mathbb{I}}-{\bf{G}})$ is actually non-invertible since it is rank-deficient; the electric potential has infinite number of solutions since it is defined up to an additive constant. This manifests from the fact that the fundamental equation we are trying to solve is the Poisson equation within the head,

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \left(\sigma {\boldsymbol{\nabla }}V\right)={\boldsymbol{\nabla }}\cdot {{\bf{J}}}^{p},\end{eqnarray} \tag{ 15 }$$

with Neumann boundary condition at each boundary (normal current continuity)

$$\begin{eqnarray}&&{\sigma }^{+}{\boldsymbol{\nabla }}V\cdot {\bf{n}}={\sigma }^{-}{\boldsymbol{\nabla }}V\cdot {\bf{n}},\end{eqnarray} \tag{ 16 }$$

where n is the outward-pointing unit surface normal. These equations are specified only up to the first derivative/gradient of V, hence V is defined up to a constant. Note that the singularity of $({\mathbb{I}}-{\bf{G}})$ is inherent to the physics of the problem at hand, and deflation is required regardless of the discretization scheme.

We thus know that both V and V + k e, where ${\bf{e}}={\left(\mathrm{1,1},\ldots ,1\right)}^{T}$ and k is a nonzero constant, are solutions to the matrix equation. So, in addition to (14), we also have

$$\begin{eqnarray}&&\left({\mathbb{I}}-{\bf{G}}\right)\left({\bf{V}}+k{\bf{e}}\right)={{\bf{V}}}_{\infty }.\end{eqnarray} \tag{ 17 }$$

Subtracting (14) from this yields

$$\begin{eqnarray}&&\left({\mathbb{I}}-{\bf{G}}\right){\bf{e}}=0\end{eqnarray} \tag{ 18 }$$

which indicates that $\left({\mathbb{I}}-{\bf{G}}\right)$ has a zero eigenvalue with associated eigenvector e ≠ 0, i.e. it is indeed singular. Equivalently,

$$\begin{eqnarray}&&{\bf{Ge}}={\bf{e}}.\end{eqnarray} \tag{ 19 }$$

This means that G has a unit eigenvalue with corresponding eigenvector e. So, one way to avoid the singularity is to eliminate this unit eigenvalue; the standard way to do this is by deflation, as follows.

First, assume the unit eigenvalue of G has a multiplicity of 1 (Lynn and Timlake 1968). For any vector a, we need to find a vector c with constant entries (not all necessarily the same) such that

$$\begin{eqnarray}&&{{\bf{c}}}^{T}{\bf{a}}=\left\{\begin{array}{l}k\ \ {\rm{if}}\ {\bf{a}}=k{\bf{e}}\\ 0\ \ {\rm{otherwise.}}\end{array}\right.\end{eqnarray} \tag{ 20 }$$

The first case imposes the condition of defining a reference potential in some way. For example, if we pick c to have all the same entries, then it means we let the sum of all potentials over all boundaries to be zero (Hämäläinen and Sarvas 1989). We may also pick just a few entries to be zero or nonzero, corresponding to a possibly more meaningful reference potential (e.g. Wilson terminals used in electrocardiogram (Fischer et al 2002)). The second case ensures that all eigenvalues of ${\bf{G}}^{\prime} =({\bf{G}}-{{\bf{ec}}}^{T})$ are equal to the eigenvalues of G, except for the unit eigenvalue which is replaced by zero. This ensures that $({\mathbb{I}}-{\bf{G}}^{\prime} )$ is non-singular and hence invertible by condition (19). We may explicitly show this preservation of eigenvalues for ${\bf{G}}^{\prime}$ as follows. Let λ and v_e be eigenvalues and corresponding eigenvectors of G. For v_e ≠ e

$$\begin{eqnarray}&&{\bf{G}}^{\prime} {{\bf{v}}}_{e}=\left({\bf{G}}-{{\bf{ec}}}^{T}\right){{\bf{v}}}_{e}={{\bf{Gv}}}_{e}-{{\bf{ec}}}^{T}{{\bf{v}}}_{e}=\lambda {{\bf{v}}}_{e}-{\bf{e}}0=\lambda {{\bf{v}}}_{e},\end{eqnarray} \tag{ 21 }$$

hence eigenvalues are preserved. For v_e = k e,

$$\begin{eqnarray}&&{\bf{G}}^{\prime} k{\bf{e}}=\left({\bf{G}}-{{\bf{ec}}}^{T}\right)k{\bf{e}}={\bf{G}}k{\bf{e}}-{{\bf{ec}}}^{T}k{\bf{e}}=k{\bf{e}}-{\bf{e}}k=0.\end{eqnarray} \tag{ 22 }$$

So, we are now able to solve for V with the invertible deflated $({\mathbb{I}}-{\bf{G}}^{\prime} )$

$$\begin{eqnarray}&&{\bf{V}}={\left({\mathbb{I}}-{\bf{G}}+{{\bf{ec}}}^{T}\right)}^{-1}{{\bf{V}}}_{\infty }.\end{eqnarray} \tag{ 23 }$$

We now have a set of discrete potential solutions at the centroids of the triangles. If we discretize the B_vol integral in an identical manner as the electric potential case above, it allows us to get an approximation of the magnetic field at arbitrary field locations r directly from (3). The magnetic field is given by

$$\begin{eqnarray}&&{\bf{B}}\left({\bf{r}}\right)\approx {{\bf{B}}}_{0}\left({\bf{r}}\right)+\displaystyle \frac{{\mu }_{0}}{4\pi }\displaystyle \sum _{l=1}^{{N}_{S}}\left({\sigma }_{l}^{-}-{\sigma }_{l}^{+}\right)\displaystyle \sum _{m=1}^{{N}_{l}}V\left({{\bf{c}}}_{l}^{m}\right){{\boldsymbol{\Omega }}}_{l}^{m}\end{eqnarray} \tag{ 24 }$$

where we have defined the "vector solid angle"

$$\begin{eqnarray}&&{{\boldsymbol{\Omega }}}_{l}^{m}\equiv {\int }_{{{\rm{\Delta }}}_{l}^{m}}\displaystyle \frac{{\bf{r}}-{\bf{r}}^{\prime} }{{\left|{\bf{r}}-{\bf{r}}^{\prime} \right|}^{3}}\times d{{\bf{S}}}_{{{\rm{\Delta }}}_{l}^{m}}^{\prime} .\end{eqnarray} \tag{ 25 }$$

The evaluation of ${{\boldsymbol{\Omega }}}_{l}^{m}$ has been done in de Munck (1992) via Stoke's Theorem, and is given by

$$\begin{eqnarray}&&{{\boldsymbol{\Omega }}}_{l}^{m}=\displaystyle \sum _{i=1}^{3}\left({\gamma }_{i-1}-{\gamma }_{i}\right){{\bf{r}}}_{i},\end{eqnarray} \tag{ 26 }$$

where

$$\begin{eqnarray}&&{\gamma }_{i}\equiv -\displaystyle \frac{1}{\left|{{\bf{r}}}_{i+1}-{{\bf{r}}}_{i}\right|}\cdot \mathrm{ln}\displaystyle \frac{{r}_{i}\left|{{\bf{r}}}_{i+1}-{{\bf{r}}}_{i}\right|+{{\bf{r}}}_{i}\cdot \left({{\bf{r}}}_{i+1}-{{\bf{r}}}_{i}\right)}{{r}_{i+1}\left|{{\bf{r}}}_{i+1}-{{\bf{r}}}_{i}\right|+{{\bf{r}}}_{i+1}\cdot \left({{\bf{r}}}_{i+1}-{{\bf{r}}}_{i}\right)},\end{eqnarray} \tag{ 27 }$$

and i = 1, 2, 3 correspond to the three vertices of triangle m on surface S_l . Also note that r₄ ≡ r₁ and r₀ ≡ r₃.

This expression allows us to calculate the magnetic field easily, assuming useful indexing of vertices and triangles has been done in the process of surface discretization.

For forward calculation of the signal vectors, we calculate the magnetic flux through magnetometer pick-up loops. In reality, this pick-up loop setup corresponds only for SQUIDs; OPMs measure the volume integral of the magnetic field over a cylindrical sensing volume. However, we are interested primarily in how the reduced distance of OPM sensor arrays may affect the signal measured due to BEM head model errors, hence we do not consider volume integrals over the magnetic field here.

For N sensors, the magnetic field is discretized into N channel readings of magnetic flux. The vectorization of these readings into an N × 1 vector ϕ is defined as the signal vector. The signal space (or signal) is the N-dimensional vector space with elements being any possible signal vector. Within the context of BEM, the signal space may be defined as the space containing all possible signals one may obtain when doing a forward calculation using the triangulation of a perfectly accurate head model.

In the case of sensors measuring the magnetic flux through a surface specified by an area and a normal vector, the jth element of ϕ is given by

$$\begin{eqnarray}&&{\phi }_{j}={\int }_{{S}_{j}}{\bf{B}}({\bf{r}})\cdot d{{\bf{S}}}_{j},\end{eqnarray} \tag{ 28 }$$

where d S_j represents an infinitesimal surface element on the sensor surface with unit normal n_j . The calculation of (28) is commonly done by cubature approximation over the sensor area (Abramowitz and Stegun 1964).

From the BEM steps above, we see that any inaccurate head mesh models will lead to inaccurate forward calculations of the magnetic flux, since they correspond to perturbed triangle vertices and centroids. This in turn leads to inaccurate calculations of the potential, magnetic field and magnetic flux signal vector. These errors may be written as an additive perturbation since they are random and independent of the unperturbed quantities,

$$\begin{eqnarray}&&{\bf{V}}^{\prime} ={\bf{V}}+\delta {\bf{V}}\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\bf{B}}^{\prime} ={\bf{B}}+\delta {\bf{B}}\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{\boldsymbol{\phi }}^{\prime} ={\boldsymbol{\phi }}+\delta {\boldsymbol{\phi }},\end{eqnarray} \tag{ 31 }$$

where the perturbations to each element of the flux signal vector are given by

$$\begin{eqnarray}&&\delta {\phi }_{j}={\int }_{{S}_{j}}\delta {\bf{B}}({\bf{r}})\cdot d{{\bf{S}}}_{j}.\end{eqnarray} \tag{ 32 }$$

Note that the BEM errors are only relevant to the forward modeling of the signal vectors; real recorded signal vectors by definition do not have BEM errors. In the context of (noiseless) BEM forward modeling, the goal is to set up a head model such that the calculated signal ${\boldsymbol{\phi }}^{\prime}$ is as close to the (noiseless) recorded/true data ϕ as possible. If information about the BEM geometry is perfect, then ${\boldsymbol{\phi }}^{\prime} ={\boldsymbol{\phi }}$ with δ ϕ = 0. We re-emphasize that the primary current contribution (4) does not depend on the head model by Geselowitz's formula, thus all errors come from inaccurate volume contribution. In other words, δ B = δ B_vol and hence δ ϕ = δ ϕ _vol .

In the Supplementary material, we present an analytical approach to calculate the first-order perturbation contributions of δ V and δ B. The subsequent calculation of the flux perturbation δ ϕ may be done numerically with cubature approximations (Abramowitz and Stegun 1964) or analytically (Yeo et al 2021).

A compact metric for quantifying the difference between the recorded/reference data and modeled/perturbed data is the angle between the corresponding signal vectors ϕ and ${\boldsymbol{\phi }}^{\prime}$ respectively. The angle between a vector y and a subspace Z can be calculated with the general formula (Gunawan et al 2005)

$$\begin{eqnarray}&&\mathrm{subspace}({\bf{p}},{\bf{Z}})=\theta =\arccos \left(\displaystyle \frac{\parallel {\mathrm{Proj}}_{{\bf{Z}}}{\bf{p}}\parallel }{\parallel {\bf{p}}\parallel }\right),\end{eqnarray} \tag{ 33 }$$

where Proj_Z p refers to the orthogonal projection of p onto Z. If the subspace Z is a vector as well, for instance, in our case where Z = ϕ and ${\bf{p}}={\boldsymbol{\phi }}^{\prime}$, then the formula reduces to the familiar $\theta =\arccos ({{\boldsymbol{\phi }}}^{T}{\boldsymbol{\phi }}^{\prime} /\parallel {\boldsymbol{\phi }}\parallel \parallel {\boldsymbol{\phi }}^{\prime} \parallel )$, and the order of the arguments in the notation 'subspace(· , · )' is interchangeable.

If we decompose the signals into their different spatial frequency components via the SSS method (Taulu and Kajola 2005), then the angle for individual different spatial frequency components can be calculated as well. In the SSS formalism, signal vectors are decomposed into their l degree components via a vector spherical harmonic (VSH) expansion, and each spatial frequency component is conveniently characterized by an l degree component. Higher l degrees correspond to higher spatial frequencies, whereas lower l degrees correspond to lower spatial frequency signal components. Let S be the basis matrix containing the VSH modes, S_1:L be the first l = L degree portion of S, and x_1:L be the corresponding coefficients/multipole moments of S_1:L . The subspace angle θ_l specific to the cumulative spatial frequency bands from l = 1 to l = L is thus

$$\begin{eqnarray}&&{\theta }_{1:L}=\arccos \left(\displaystyle \frac{\parallel {\mathrm{Proj}}_{{{\boldsymbol{\phi }}}_{1:L}}{{\boldsymbol{\phi }}}_{1:L}^{\prime} \parallel }{\parallel {{\boldsymbol{\phi }}}_{1:L}^{{\prime} }\parallel }\right),\end{eqnarray} \tag{ 34 }$$

where

$$\begin{eqnarray}&&{{\boldsymbol{\phi }}}_{1:L}={{\bf{S}}}_{1:L}{{\bf{x}}}_{1:L}\end{eqnarray} \tag{ 35 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{\phi }}}_{1:L}^{{\prime} }={{\bf{S}}}_{1:L}{{\bf{x}}}_{1:L}^{{\prime} }.\end{eqnarray} \tag{ 36 }$$

Given a measured signal vector ϕ , the multipole moments x_1:L must first be estimated by taking the pseudoinverse of S_1:L ,

$$\begin{eqnarray}&&{{\bf{x}}}_{1:L}\approx {{\bf{S}}}_{1:L}^{\dagger }{\boldsymbol{\phi }}.\end{eqnarray} \tag{ 37 }$$

Then, if required, individual l degree portions x_1:L may be obtained from x_1:L . Note that to avoid aliasing in signal reconstruction, a high enough l degree truncation is required so that the high frequency components of the signal do not get projected inaccurately onto basis vectors corresponding to low frequencies. It was shown in Taulu and Kajola (2005) that a truncation at approximately L = 8 is sufficient to represent signal vectors; in this paper, we truncate at L = 12 in anticipation of close sensor array distances capturing signals of higher spatial complexities.

For our reference set-up, we used a simple 1-shell spherical head model of radius 9 cm and conductivity 0.33 S m⁻¹ with origin located at the center of the sphere, and did a CC BEM forward calculation as described in sections 2 and 3 to obtain ϕ . The BEM method was implemented using the Matlab library as provided in Stenroos et al (2007). The spherical surface was triangulated with 642 vertices, and 4 dipole sources with varying distances located at (2, 0, 0) cm, (4, 0, 0) cm, (6, 0, 0) cm, and (8, 0, 0) cm from the origin were specified. All the dipoles had a moment of (0, 10, 0) nAm. To simulate inaccurate mesh modeling to obtain ${\boldsymbol{\phi }}^{\prime}$, the vertices of the spherical mesh were randomly perturbed radially up to 2%, 4%, 6%, 8% and 10% relative to the 9 cm radius.

The sensor array that we used consisted of 324 square magnetometer pick-up loops with side length 2.1 cm all with non-radial orientations (to avoid linear dependence of the SSS basis (Taulu and Kajola 2005)), uniformly arranged on a spherical shell up to π/6 below the z = 0 plane. It has been shown in Yeo et al (2021) that a 9-point cubature approximation yields accurate calculations for the sensor distances of 10 cm to 15 cm that we considered, thus we used a 9-point cubature approximation for the calculation of signal vectors in this paper. Figure 1 illustrates the sensor, head, and source setup as described above.

Zoom In Zoom Out Reset image size

The averaged results over 100 forward calculations with this random vertex perturbation setup is presented in the following sections.

First, we investigate how the error due to BEM mesh inaccuracy varies according to sensor array distance using equation (33). The sensor array radii were set to be from 10 cm to 15 cm, in increments of 1 cm (i.e. 1 cm to 6 cm from the surface of a 9 cm head model) and the signal was assumed the be noiseless. Figure 2 shows that for all the source distances considered, as sensor array distance increases, the subspace angle between the reference signal ϕ and perturbed signal ${\boldsymbol{\phi }}^{\prime}$ decreases, indicating decreasing relative effects of mesh boundary inaccuracies as sensor array distance from the head increases. Moreover, smaller perturbations to the head model resulted in smaller subspace angles for a given sensor array distance when compared to higher perturbations as expected. These results may also be visually seen via plots of ϕ and ${\boldsymbol{\phi }}^{\prime}$ explicitly; we show this in figure 3 with the 2 cm source case across the various sensor distances and mesh perturbations for one of the 100 random mesh perturbation realizations. Also note the interesting observation from figure 2 that there is a 'turnaround' point in subspace angle as a function of source distance for a fixed mesh perturbation amount—from the 2 cm source case to the 4 cm and 6 cm source cases, the subspace angle decreases before increasing again in the 8 cm source case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we determine if the additive error of the signal vector, δ ϕ , may be explained by increasing orders of l degree truncation, and if δ ϕ varies according to sensor distances. First, we show that the unperturbed signal may be explained to a large extent with L = 12 degree truncation for all the sensor array distances that we considered. This is shown in the first row of plots in figure 4; the subspace angle with S_1:L decays to become small (<5°) at L = 12 for all source distances. As expected, closer sensor array distances have higher subspace angles than further sensor arrays, which agrees with figure 2.

Zoom In Zoom Out Reset image size

From the second row of plots in figure 4 however, we see that despite the total signal ϕ being well-explained at L = 12, the signal errors δ ϕ which consist mostly of higher frequency components require a higher L truncation to explain the signal for all source distances. The subspace angle still decreases for increased l degree truncation as expected (since δ ϕ = δ ϕ_vol are still volume current contributions), and the subspace angles are smaller for increased sensor array distances. However, the angles are much higher than those of the total signal ϕ , indicating that higher spatial frequency components are affected more by inaccurate mesh modeling. This is because a higher l degree truncation required to fully explain the additive signal error δ ϕ of closer sensor arrays indicates that their δ ϕ has a more spatially complex pattern with higher spatial frequency components. This higher spatial complexity is due to the fact that higher spatial frequency components decay quickly with increasing distance, hence by placing sensors closer to the source, these components can be captured to give higher resolution of the signal error δ ϕ . Note that figure 4 is for the case where the vertices of the mesh was perturbed radially up to 10%. The plots for the different mesh perturbations are similar, hence we present just the 10% case here.

From our results above which considers the noiseless signal case, the signal vector suffers from higher inaccuracies for closer sensor arrays due to the larger errors in their higher spatial frequency components. Deep and superficial sources also had higher signal inaccuracies due to a higher volume current contribution relative to primary currents. Here, we investigate if these observations will be seen in the form of source localization errors as well. The source localization procedure was done via a standard ECD fit using the 'fit_dipole' function in MNE-Python 1.0 (Gramfort et al 2013, 2014). The forward model here was computed using LC BEM however, as it is currently the only solving method implemented in MNE-Python.

The results for the noiseless case are shown in the first columns of figures 5, 6, 7 and 8, which correspond to the 2 cm, 4 cm, 6 cm and 8 cm source case respectively. Indeed, in the noiseless case, source localization and orientation errors are largest for closer sensor arrays, as well as deep and superficial sources for a fixed mesh perturbation. These observations are in agreement with the result for signal errors as presented above.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, the primary motivation to move sensors closer to the head is because closer sensors have the potential to give more accurate source localization results due to higher SNR as mentioned in the introduction. The definition of SNR (with units of decibels) used in MNE-Python 1.0 follows that of Goldenholz (Goldenholz et al 2008)

$$\begin{eqnarray}&&\mathrm{SNR}=10\mathrm{log}\unicode{x0230A}\displaystyle \frac{{q}^{2}}{N}\displaystyle \sum _{j=1}^{N}\displaystyle \frac{{\phi }_{j}^{2}}{{s}_{j}^{2}}\unicode{x0230B},\end{eqnarray} \tag{ 38 }$$

where q is the source amplitude, N is the total number of sensors in the sensor array, ϕ_j is the signal on sensor j for a source with unit amplitude, and s_j ² is the noise variance on sensor j. We thus introduced various noise levels to determine the noise level at which the effect of having more accurate source localization due to higher SNR of closer sensors outweighs the effect of less accurate source localization due to the higher errors of the high frequency components of the noiseless signal. We found that this occurred at a SNR of around 6 dB (the sensor noise was varied to force a 6 dB SNR for a varying sensor array distances); the localization and orientation errors began to increase as sensor array distances increased at around an SNR of 6 dB. This is shown in the second column of figures 5, 6, 7 and 8. The rightmost column of figures 5, 6, 7 and 8 show that for an SNR greater than 6 dB, in this case a constant 20 fT noise level, then there is an improvement in source localization for closer sensor arrays. This means that a higher SNR allowing for better source localization resolutions outweighs the effects of BEM errors for noisy signals with SNR above approximately 6 dB.