Integral equations and boundary-element solution for static potential in a general piece-wise homogeneous volume conductor

Study a compartment V_i of conductivity ${{\sigma}_{i}}$ bounded by closed surface $\partial {{V}_{i}}$. Label compartments and surfaces using sub- and superscripts, respectively. Using Green's theorem and (1), we get (Stenroos 2009)

$$\begin{eqnarray}&&{{ \Omega }^{\partial {{V}_{i}}}}\left(\boldsymbol{r}\right)\phi \left(\boldsymbol{r}\right)={{v}_{i}}\left(\boldsymbol{r}\right)+\left({{S}^{\partial {{V}_{i}}}} \Gamma _{-}^{\partial {{V}_{i}}}\right)\left(\boldsymbol{r}\right)-\left({{D}^{\partial {{V}_{i}}}}\phi \right)\left(\boldsymbol{r}\right),\end{eqnarray} \tag{ 3 }$$

where ${{ \Omega }^{\partial {{V}_{i}}}}\left(\boldsymbol{r}\right)$ is the normalized solid angle spanned by the boundary of V_i at $\boldsymbol{r}$, v_i is the potential produced in an infinite homogeneous medium of conductivity ${{\sigma}_{i}}$ by sources within V_i, Sⁱ is the single-layer operator, and D^j is the double-layer operator:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{v}_{i}}\left(\boldsymbol{r}\right)=\frac{1}{4\pi {{\sigma}_{i}}}{{{\int}^{}}_{{{V}_{i}}}}\frac{{{i}_{\text{v}}}\left(\boldsymbol{r}^{{\prime}}\right)}{|\boldsymbol{r}-\boldsymbol{r}^{{\prime}}|}\,\text{d}{{V}^{\prime}} \end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \left({{D}^{j}}f\right)\left(\boldsymbol{r}\right)=\frac{1}{4\pi}{{{\int}^{}}_{{{S}^{j}}}}f\left(\boldsymbol{r}^{{\prime}}\right)\frac{\left(\boldsymbol{r}-\boldsymbol{r}^{{\prime}}\right)}{|\boldsymbol{r}-\boldsymbol{r}^{{\prime}}{{|}^{3}}}\centerdot \,{{\overrightarrow{\text{d}S}}^{\prime}} \end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \left({{S}^{j}}f\right)\left(\boldsymbol{r}\right)=\frac{1}{4\pi}{{{\int}^{}}_{{{S}^{j}}}}\frac{f\left(\boldsymbol{r}^{{\prime}}\right)}{|\boldsymbol{r}-\boldsymbol{r}^{{\prime}}|}\,\text{d}{{S}^{\prime}}. \end{array}\end{eqnarray} \tag{ 6 }$$

Now consider a two-compartment (${{V}_{a}},{{V}_{b}}$) volume conductor immersed in infinite homogeneous space V_c as illustrated figure 2. Compartments a and b are connected with each other via surface S¹ and with compartment c via surfaces S² and S³, respectively. All surfaces are open and join each other at junctions, which form the boundaries of the surfaces. Applying (3) to compartments a,b and c, we get

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{ \Omega }^{\partial {{V}_{a}}}}\phi \left(\boldsymbol{r}\right)={{v}_{a}}\left(\boldsymbol{r}\right)+\left({{S}^{1}} \Gamma _{-}^{1}\right)\left(\boldsymbol{r}\right)+\left({{S}^{2}} \Gamma _{-}^{2}\right)\left(\boldsymbol{r}\right)-\left({{D}^{1}}{{\phi}^{1}}\right)\left(\boldsymbol{r}\right)-\left({{D}^{2}}{{\phi}^{2}}\right)\left(\boldsymbol{r}\right) \end{array}\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{ \Omega }^{\partial {{V}_{b}}}}\phi \left(\boldsymbol{r}\right)={{v}_{b}}\left(\boldsymbol{r}\right)-\left({{S}^{1}} \Gamma _{+}^{1}\right)\left(\boldsymbol{r}\right)+\left({{S}^{3}} \Gamma _{-}^{3}\right)\left(\boldsymbol{r}\right)+\left({{D}^{1}}{{\phi}^{1}}\right)\left(\boldsymbol{r}\right)-\left({{D}^{3}}{{\phi}^{3}}\right)\left(\boldsymbol{r}\right) \end{array}\end{eqnarray} \tag{ 7b }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{ \Omega }^{\partial {{V}_{c}}}}\phi \left(\boldsymbol{r}\right)={{v}_{c}}\left(\boldsymbol{r}\right)-\left({{S}^{2}} \Gamma _{+}^{2}\right)\left(\boldsymbol{r}\right)-\left({{S}^{3}} \Gamma _{+}^{3}\right)\left(\boldsymbol{r}\right)+\left({{D}^{2}}{{\phi}^{2}}\right)\left(\boldsymbol{r}\right)+\left({{D}^{3}}{{\phi}^{3}}\right)\left(\boldsymbol{r}\right). \end{array}\end{eqnarray} \tag{ 7c }$$

Zoom In Zoom Out Reset image size

Next, multiply each equation by the conductivity of the corresponding compartment. For surface terms, express the conductivity in terms of the corresponding surface, labeling compartments inside and outside the surface k as $\sigma _{-}^{k}$ and $\sigma _{+}^{k}$, respectively; for example, ${{\sigma}_{b}}=\sigma _{+}^{1}=\sigma _{-}^{3}$. Sum the equations and use (2b) to get

$$\begin{eqnarray}&&\underset{i=a}{\overset{c}{\sum}}\,{{\sigma}_{i}}{{ \Omega }^{\partial {{V}_{i}}}}\phi \left(\boldsymbol{r}\right)=\underset{i=a}{\overset{c}{\sum}}\,{{\sigma}_{i}}{{v}_{i}}\left(\boldsymbol{r}\right)+\underset{l=1}{\overset{3}{\sum}}\,\left(\sigma _{+}^{l}-\sigma _{-}^{l}\right)\left({{D}^{l}}{{\phi}^{l}}\right)\left(\boldsymbol{r}\right).\end{eqnarray} \tag{ 8 }$$

The first right-side term results in a conductivity-independent potential function; we express it using infinite-medium potential ${{\phi}_{\infty}}$ in a medium that has dummy (unit) conductivity ${{\sigma}_{\text{s}}}$. Equation (8) generalizes directly to M compartments and N boundary surfaces, resulting in the general form of surface integral equation for scalar potential:

$$\begin{eqnarray}&&\underset{i=1}{\overset{M}{\sum}}\,{{\sigma}_{i}}{{ \Omega }^{\partial {{V}_{i}}}}\phi \left(\boldsymbol{r}\right)={{\sigma}_{\text{s}}}{{\phi}_{\infty}}\left(\boldsymbol{r}\right)+\frac{1}{4\pi}\underset{l=1}{\overset{N}{\sum}}\,\left(\sigma _{+}^{l}-\sigma _{-}^{l}\right){\int}_{{{S}^{l}}}\phi \left(\boldsymbol{r}^{{\prime}}\right)\frac{\left(\boldsymbol{r}-\boldsymbol{r}^{{\prime}}\right)}{|\boldsymbol{r}-\boldsymbol{r}^{{\prime}}{{|}^{3}}}\centerdot \,{\vec{\rm d}}{{S}^{\prime}},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{{\phi}_{\infty}}\left(\boldsymbol{r}\right)=\frac{1}{4\pi {{\sigma}_{\text{s}}}}{{{\int}^{}}_{\text{all}~{{i}_{\text{v}}}}}\frac{{{i}_{\text{v}}}\left(\boldsymbol{r}^{{\prime}}\right)}{|\boldsymbol{r}-\boldsymbol{r}^{{\prime}}|}\,\text{d}{{V}^{\prime}}=\frac{1}{4\pi {{\sigma}_{\text{s}}}}{\int}_{\text{all}~{{{\vec{J}}}_{\text{s}}}}\frac{{{{\vec{J}}}_{\text{s}}}\left(\boldsymbol{r}^{{\prime}}\right)\centerdot \left(\boldsymbol{r}-\boldsymbol{r}^{{\prime}}\right)}{|\boldsymbol{r}-\boldsymbol{r}^{{\prime}}|}\,\text{d}{{V}^{\prime}}.\end{eqnarray} \tag{ 10 }$$

Evaluating the left side of (9) with different positionings of $\boldsymbol{r}$ but avoiding junctions, the equation reduces to previously presented equations.

• For $\boldsymbol{r}\in {{V}_{j}},\boldsymbol{r}\notin \partial {{V}_{j}}$, we have ${{ \Omega }^{\partial {{V}_{j}}}}\left(\boldsymbol{r}\right)=1$ and ${{ \Omega }^{\partial {{V}_{k}}}}\left(\boldsymbol{r}\right)=0$, $j\ne k$, yielding the equation presented by Geselowitz (1967).
• For $\boldsymbol{r}\in {{S}^{k}}$ and $\boldsymbol{r}\notin {{S}^{j}},j\ne k$, the left side yields $\left(\sigma _{+}^{k}{{ \Omega }^{\partial V_{+}^{k}}}+\sigma _{-}^{k}{{ \Omega }^{\partial V_{-}^{k}}}\right)\phi \left(\boldsymbol{r}\right)$, where $V_{\pm }^{k}$ label compartments outside and inside S^k, respectively. The result matches the sharp-edged form (Ferguson and Stroink 1997, Stenroos 2009).
• If $\boldsymbol{r}\in {{S}^{k}}$ and $\boldsymbol{r}\notin {{S}^{j}},j\ne k$, and S^k is smooth around $\boldsymbol{r}$, we get the form presented by Barnard et al (1967b) and in junctioned geometry by Hämäläinen et al (1993),
  $$\begin{eqnarray}&&\frac{1}{2}\left(\sigma _{+}^{k}+\sigma _{-}^{k}\right)\phi \left(\boldsymbol{r}\right)={{\sigma}_{\text{s}}}{{\phi}_{\infty}}\left(\boldsymbol{r}\right)+\frac{1}{4\pi}\underset{l=1}{\overset{N}{\sum}}\,\left(\sigma _{+}^{l}-\sigma _{-}^{l}\right){\int}_{{{S}^{l}}}\phi \left(\boldsymbol{r}^{{\prime}}\right)\frac{\left(\boldsymbol{r}-\boldsymbol{r}^{{\prime}}\right)}{|\boldsymbol{r}-\boldsymbol{r}^{{\prime}}{{|}^{3}}}\centerdot \,{\vec{\rm d}}{{S}^{\prime}}.\end{eqnarray} \tag{ 11 }$$

A boundary-element model is built by tessellating boundary surfaces into polygon meshes, approximating potential ϕ on all boundary surfaces as a linear combination of a set of basis functions defined around a total of N_d discretization points, evaluating (9) on all boundary surfaces, and minimizing the residual of the approximated potential with respect to N_d weight functions. This results in an equation of the form $\boldsymbol{T\Phi}=\boldsymbol{}{{\mathbf{\Phi}}^{\infty}}$, where $\left({{N}_{\text{d}}}\times 1\right)$-vectors $\boldsymbol{\Phi}$ and $\boldsymbol{}{{\mathbf{\Phi}}^{\infty}}$ contain the values of ϕ and weighted ${{\sigma}_{\text{s}}}{{\phi}_{\infty}}$ in all discretization points, and $\mathbf{T}$ is a $\left({{N}_{\text{d}}}\times {{N}_{\text{d}}}\right)$ matrix. The process is thoroughly explained in, e.g. Stenroos et al (2007), Stenroos (2009), and the assembly of $\mathbf{T}$ and $\boldsymbol{}{{\mathbf{\Phi}}^{\infty}}$ as used in this work is shown compactly in the appendix.

In this work, I use triangle meshes, linear basis functions, and collocation and Galerkin weightings (de Munck 1992, Mosher et al 1999): triangle mesh k that contains $N_{\text{v}}^{k}$ vertices and $N_{\text{t}}^{k}$ triangles is defined by an $\left(N_{\text{v}}^{k}\times 3\right)$ array of vertex coordinates and an $\left(N_{\text{t}}^{k}\times 3\right)$ array of indices that point to rows of the vertex array ('local indexing'), and a linear basis function is spanned on all triangles that belong to one vertex $\boldsymbol{v}_{i}^{k}$ so that the function gets value 1 in the target vertex and decreases linearly to value 0 in the neighboring vertices. N_d is thus equal to the total number of vertices ${\sum}_{k=1}^{N}N_{\text{v}}^{k}$ in the model. In collocation weighting, Dirac delta functions defined in vertices are used as weight functions, while in Galerkin weighting, the same functions are used as basis and weight functions (Mosher et al 1999, Stenroos and Haueisen 2008). In $\boldsymbol{\Phi},\boldsymbol{}{{\mathbf{\Phi}}^{\infty}}$ and $\mathbf{T}$, vertices of each mesh are concatenated together so that the vertices of kth mesh have pooled indices of $\left[\left({\sum}_{i=1}^{(k-1)}N_{\text{v}}^{i}\right)+1,{\sum}_{i=1}^{k}N_{\text{v}}^{i}\right]$.

In non-junctioned geometry, meshes are closed and fully separate; every row of $\boldsymbol{\Phi}$ and $\boldsymbol{}{{\mathbf{\Phi}}^{\infty}}$ refers to a unique vertex, and every element of $\mathbf{T}$ refers to a unique pair of two vertices. In junctioned geometry, the meshes that are connected to a junction are open, and vertices and element edges at junctions are shared by at least two meshes. Using standard BEM implementation with these meshes, many rows of $\boldsymbol{\Phi}$ and $\boldsymbol{}{{\mathbf{\Phi}}^{\infty}}$ and many elements of $\mathbf{T}$ thus refer to the same vertices or vertex pairs. Now consider a junctioned geometry. The standard BEM matrix is of the form

$$\begin{eqnarray}&&{{\mathbf{T}}_{\text{p}}}\boldsymbol{}{{\mathbf{\Phi}}_{\text{p}}}=\boldsymbol{\Phi}_{\text{p}}^{\infty},\end{eqnarray} \tag{ 12 }$$

where subscript p refers to pooled indexing. Collect all unique vertices to a global set of N_u vertices, ${{N}_{\text{u}}}<{{N}_{\text{d}}}$. Then define matrix operators that convert from pooled to global indexing in either summing or extracting way: operator ${{\mathbf{I}}_{\text{p}2\text{g}}}$ (${{N}_{\text{u}}}\times {{N}_{\text{d}}}$ matrix) converts data in rows from pooled to global indexing and sums the contributions from junction vertices, while ${{\mathbf{E}}_{\text{p}2\text{g}}}$ does the same but extracts the contribution from only one of the junction vertices. Operator ${{\mathbf{I}}_{\text{g}2\text{p}}}$ converts from global to pooled indexing, setting the global values to all corresponding vertices. It turns out that ${{\mathbf{I}}_{\text{g}2\text{p}}}=\mathbf{I}_{\text{p}2\text{g}}^{\text{T}}$, where ^T marks transposis.

In linear collocation (LC) BEM, each row of (12) refers to one vertex. The rows that map to the same global vertex contain the same information, and we can use the ${{\mathbf{E}}_{\text{p}2\text{g}}}$ operator to convert to global indexing:

$$\begin{eqnarray}&&{{\mathbf{E}}_{\text{p}2\text{g}}}{{\mathbf{T}}_{\text{p}}}\boldsymbol{}{{\mathbf{\Phi}}_{\text{p}}}={{\mathbf{E}}_{\text{p}2\text{g}}}\boldsymbol{\Phi}_{\text{p}}^{\infty}=\boldsymbol{\Phi}_{\text{g}}^{\infty}.\end{eqnarray} \tag{ 13 }$$

Each column of ${{\mathbf{T}}_{\text{p}}}$ models the contribution of the neighborhood of a vertex. For junction vertices, these neighborhoods are different for different meshes. To convert ${{\mathbf{T}}_{\text{p}}}$ to operate on global indexing, we sum the contributions of each local vertex using ${{\mathbf{I}}_{\text{p}2\text{g}}}$. To apply ${{\mathbf{I}}_{\text{p}2\text{g}}}$ to columns, transpose it and multiply ${{\mathbf{T}}_{\text{p}}}$ from the right side,

$$\begin{eqnarray}&&{{\mathbf{E}}_{\text{p}2\text{g}}}{{\mathbf{T}}_{\text{p}}}\mathbf{I}_{\text{p}2\text{g}}^{\text{T}}\boldsymbol{}{{\mathbf{\Phi}}_{\text{g}}}=\boldsymbol{\Phi}_{\text{g}}^{\infty}.\end{eqnarray} \tag{ 14 }$$

Another way to explain this step is to just say that we express $\boldsymbol{}{{\mathbf{\Phi}}_{\text{p}}}$ as ${{\mathbf{I}}_{\text{g}2\text{p}}}\boldsymbol{}{{\mathbf{\Phi}}_{\text{g}}}=\mathbf{I}_{\text{p}2\text{g}}^{\text{T}}\boldsymbol{}{{\mathbf{\Phi}}_{\text{g}}}$. We now get the final transfer matrix ${{\mathbf{F}}_{\text{g}}}$ in global indexing:

$$\begin{eqnarray}&&\boldsymbol{}{{\mathbf{\Phi}}_{\text{g}}}=\underset{{{\mathbf{F}}_{\text{g}}}}{{\underbrace{{{{\left({{\mathbf{E}}_{\text{p}2\text{g}}}{{\mathbf{T}}_{\text{p}}}{{\mathbf{I}}_{\text{g}2\text{p}}}\right)}^{-1}}}}\,}}\,\boldsymbol{\Phi}_{\text{g}}^{\infty}.\end{eqnarray} \tag{ 15 }$$

For linear Galerkin (LG) BEM, the logic is the same, but as also rows characterize neighborhoods, we need to use the summing index conversion:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathbf{I}}_{\text{p}2\text{g}}}{{\mathbf{T}}_{\text{p}}}\mathbf{I}_{\text{p}2\text{g}}^{\text{T}}\boldsymbol{}{{\mathbf{\Phi}}_{\text{g}}}={{\mathbf{I}}_{\text{p}2\text{g}}}\boldsymbol{\Phi}_{\text{p}}^{\infty} \end{array}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Rightarrow \boldsymbol{}{{\mathbf{\Phi}}_{\text{g}}}=\underset{{{\mathbf{F}}_{\text{p}2\text{g}}}}{{\underbrace{{{{\left({{\mathbf{I}}_{\text{p}2\text{g}}}{{\mathbf{T}}_{\text{p}}}\mathbf{I}_{\text{p}2\text{g}}^{\text{T}}\right)}^{-1}}{{\mathbf{I}}_{\text{p}2\text{g}}}}}\,}}\,\boldsymbol{\Phi}_{\text{p}}^{\infty}. \end{array}\end{eqnarray} \tag{ 17 }$$

Here $\boldsymbol{\Phi}_{\text{p}}^{\infty}$ is kept because it can be directly implemented using standard BEM routines and local meshes, while implementing the LG version of $\boldsymbol{\Phi}_{\text{g}}^{\infty}$ would demand the construction of a junctioned global mesh, adding unnecessary complexity.

I implemented the index conversion operators and formulas (15) and (17) in Matlab (version R2014b, www.mathworks.com) using existing BEM tools (Stenroos et al 2007, Stenroos and Haueisen 2008, Stenroos 2009) and carried out three types of verifications:

First, I verified the index conversions by randomly splitting the boundary of a spherical single-compartment model to three open meshes, building BEM models using both the original closed mesh (model 1) and the separate open meshes (model 2), and simulating a set of random dipolar test sources; the resulting surface potentials in models 1 and 2 were identical up to numerical precision (relative difference RDIF¹  ∼10⁻¹⁵). Next, I split a spherical model to two open half-spheres at xy plane and added a boundary between the halfs, resulting in a two-compartment three-mesh model as illustrated in figure 3. First, I set ${{\sigma}_{1}}={{\sigma}_{2}}$ (model 3), resulting in a model physically identical to model 1, and evaluated the outer-boundary potential of a random set of test sources. With the LC approach, the results of model 3 should be nearly identical to model 1, and indeed they were (RDIF  ∼10⁻¹⁵). With LG approach, differences are expected to occur near the junctioned boundary at z  =  0, because the potential for the junction vertices is integrated in different neighborhoods in models 1 and model 3. These differences were, however, much smaller than overall errors of these models as compared to the closed-form analytical formula.

Zoom In Zoom Out Reset image size

To verify the weighting of conductivity terms, I set ${{\sigma}_{1}}$ to 1 S m⁻¹ and ${{\sigma}_{2}}$ to 0.5 S m⁻¹ or 0.1 S m⁻¹, and implemented the same models using Comsol Multiphysics software (version 5.1, www.comsol.com) that uses the finite element method (FEM). The simulation scenario had to be easy for the FEM to handle and clearly show the effect of the conductivity boundary. To achieve this, I placed a monopolar sink and source symmetrically around the origin, rather close to the split but far from each other (source locations  ±(8, 8, 2) cm, strengths  ±0.1 A, sphere radius 13 cm), and evaluated potential on the xz plane; the BEM solution for the volume was obtained by first solving potential on all boundary surfaces and then using (9) with field points set in the evaluation plane. The results are shown in figure 3. The field patterns obtained using the highest-resolution physics-controlled mesh in Comsol ('extremely fine', degrees of freedom $\text{DoF}=1542\,049$) and BEM ($\text{DoF}=3205$) were nearly identical ($0.001<\text{RDIF}<0.003$), suggesting that the BEM solver is correctly implemented. The difference between Comsol 'extremely fine' and 'fine' ($\text{DoF}=23\,061$) meshings was slightly larger ($0.009<\text{RDIF}<0.018$).

Consider a four-shell (brain, cerebrospinal fluid CSF, skull, scalp) volume conductor, where the skull has a small hole that is filled with CSF. In this simulation, I compare surface potentials obtained using FEM and four Galerkin-based BEM approaches. I built a spherical concentric four-shell model (radii 78, 81, 85, and 89 mm, conductivities 0.33, 1.79, 0.0132, and 0.33 S m⁻¹), manually made a cavity (mean diameter 23 mm, opening angle  ≈${{16}^{\circ}}$) through the skull compartment, and included the cavity in the CSF compartment. The model geometry is illustrated figure 4(a). I implemented the model in Comsol using default settings for meshing and solver (normal physics-controlled meshing, quadratic solver (QFEM, $\text{DoF}=102\,802$) and also a linear solver (LFEM, $\text{DoF}=13\,470$) and extracted the boundary meshes for the BEM models (number of vertices $=\text{DoF}=9801$). The LFEM had the same surface elements and approximation order as the Galerkin BEM and thus enables direct and fair comparison of surface and volume approaches. Using a small set of test dipoles, I assessed the error of the QFEM as function of source depth against a denser FEM model ($\text{DoF}=240\,329$) and considered QFEM an adequate reference model, when sources are at max. radius of 72 mm, where the median RDIF between QFEM and the denser FEM was 1%.

Zoom In Zoom Out Reset image size

I generated a random set of 100 source positions on a spherical cap that was centered around the cavity and had the opening angle of 32°; approximately one third of the sources were under the cavity. For these template sources, I assigned randomly-oriented unit-strength dipole moments and projected the sources at 10 different depths, totaling in 1000 test sources. I solved the surface potentials of all sources using the QFEM, LFEM and BEMs built with four different linear approaches: the standard Galerkin (LG) with 13 quadrature points per triangle, the simpler and faster quick Galerkin (LGQ) (Stenroos and Nenonen 2012) that has only one quadrature point per triangle, and isolated-source formulations of these (LGISA, LGQISA) that used approach 3 of Stenroos and Sarvas (2012). Then, I computed relative differences between the QFEM model and the other models. The results as function of source depth are shown in figure 4(b). LG and LGISA perform identically well, with median RDIF below 1.2% for all source depths. LGQISA stays below 1% at depths below 60 mm; with more superficial sources the error increases to about 2.2%. LGQ without ISA is clearly less accurate, when sources approach the pial boundary and skull. LFEM is the least accurate of the tested approaches. The maximum RDIFs of BEMs and LFEM were approximately 1.5–2.2 and 2.8 times the median RDIF, respectively.

Finally, I illustrate the effect of such a hole on electroencephalographic scalp potentials: I placed tangential dipoles (dipole moment 100 nAm) at the depth of 10 mm from the CSF boundary at three positions (1) far from the hole at 45° angle from the center of the hole; (2) under the hole boundary at 8°, and (3) under the hole center at 0°. The scalp potentials are shown in figures 4(c)–(e). For the source at 45° (c), the scalp topography is highly similar to that produced without the hole. When the source is under the boundary of the hole (d), the hole strongly affects the shape of the topography and the overall effect of the hole is at largest.