Sensitivity of the Projected Subtraction Approach to Mesh Degeneracies and Its Impact on the Forward Problem in EEG

Objective: Subtraction-based techniques are known for being theoretically rigorous and accurate methods for solving the forward problem in electroencephalography (EEG-FP) by means of the finite-element method. Within them, the projected subtraction (PS) approach is generally adopted because of its computational efficiency. Although this technique received the attention of the community, its sensitivity to degenerated elements is still poorly understood. In this paper, we investigate the impact of low-quality tetrahedra on the results computed with the PS approach. Methods: We derived upper bounds on the relative error of the element source vector as a function of geometrical features describing the tetrahedral discretization of the domain. These error bounds were then utilized for showing the instability of the PS method with regards to the mesh quality. To overcome this issue, we proposed an alternative technique, coined projected gradient subtraction (PGS) approach, that exploits the stability of the corresponding bounds. Results: Computer simulations showed that the PS method is extremely sensitive to the mesh shape and size, leading to unacceptable solutions of the EEG-FP in case of using suboptimal tessellations. This was not the case of the PGS approach, which led to stable and accurate results in a comparable amount of time. Conclusion: Solutions of the EEG-FP computed with the PS method are highly sensitive to degenerated elements. Such errors can be mitigated by the PGS approach, which showed better performance than the PS technique. Significance: The PGS is an efficient method for computing high-quality lead field matrices even in the presence of degenerated elements.

Thi s v e r sio n is b ei n g m a d e a v ail a bl e in a c c o r d a n c e wit h p u blis h e r p olici e s. S e e h t t p://o r c a . cf. a c. u k/ p olici e s. h t ml fo r u s a g e p olici e s. Co py ri g h t a n d m o r al ri g h t s fo r p u blic a tio n s m a d e a v ail a bl e in ORCA a r e r e t ai n e d by t h e c o py ri g h t h ol d e r s .
problems, and its achievement is often considered as a guarantee of accurate outcomes. However, only suboptimal elements are reachable in realistic situations comprising convoluted geometries. For this reason, gaining insights into the relation between the mesh geometry and the numerical accuracy of a technique is fundamental for understanding its sensitivity to the mesh, and consequently a key indicator of its robustness [1]. In this paper we present a detailed analysis on the relation between the mesh geometry and the numerical accuracy in the solution of the forward problem in electroencephalography (EEG-FP). This problem consists in a Poisson-like equation (subject to a Neumann boundary condition) representing the electric potential distribution in the head generated by a set of known sources of brain activity [2]. These sources are generally modelled as dipoles, introducing a singularity into the differential formulation. Within all the existing FE techniques dealing with such singularities, we focus our attention in the subtraction approach [3], [4]. This methodology stands out of the rest for allowing the use of truly (i.e., non-approximated) dipolar (or either multipolar) sources while guaranteeing the existence and uniqueness of the solution. Subtraction-based techniques were shown to provide highly accurate solutions to the EEG-FP, improving those obtained by the partial integration method, and comparable to using the Saint Venant's principle [5]- [7]. The major constraint found by this approach is related to the prohibitive computational efforts required for its use. This led Wolters et al. [3] to present the projected subtraction (PS) approach, in which the computational cost is highly reduced at the expense of accuracy. Then, it is of interest to generate accurate and fast approximations to translate the theoretical advantages of subtraction techniques to the practice. Moreover, previous work on subtraction methods has only focused on the characterisation of discretisation errors due to approximations in the source vector. However, the impact of the mesh quality in the corresponding solutions is still poorly understood.
The analysis presented here is based on the use of tetrahedral meshes and linear FE basis functions as this is the most common and widely accepted procedure in the field (e.g., [3], [4], [8]). In this context, we show that the PS approach is extremely sensitive to the mesh geometry, and therefore not recommended when the tetrahedral discretisation is not of highest quality. The reason for this is that the PS reduces its computational load by approximating the gradient of the non-linear singularity function over the mesh with the gradient of the interpolated (linear) function, whose error bounds are known to be extremely sensitive to the elements' shape [1]. To overcome this issue, we propose a new approach in which the gradient of the singularity function is approximated by the (linear) interpolation of the gradients on the nodes, exploiting the stable nature of the corresponding bounds [1]. The resulting method, coined projected gradient subtraction (PGS) approach, allows to mitigate the sensitivity of the solution to the quality of the tetrahedral mesh inherent to the PS, while preserving its computational efficiency. We demonstrate that such approximation turns the method more stable and robust to the presence of degenerated tetrahedra, a more than frequent situation when working with mesh generators based on the Delaunay triangulation [4], [9]- [11]. The applicability of the method is illustrated in the solution of the inverse problem in EEG (EEG-IP), whose accuracy relies on the technique utilised to solve the EEG-FP.

1) Differential Formulation:
The EEG-FP consists in finding the electric potential function u(r) due to a current source with density s(r) defined over the domain Ω (i.e., the head), with boundary Γ.L e tσ(r) be the rank-2 conductivity tensor field within Ω, andň(r) the unitary vector normal to Γ. Then, under generally accepted assumptions (as the quasistatic and the point electrode model approximations), the EEG-FP reduces to find u(r) satisfying ∇· σ(r)∇u(r) = −s(r) (r ∈ Ω), together with the boundary condition σ(r)∇u(r),ň(r) =0 (r ∈ Γ) [2]. In the case of assuming a dipolar source located in r 0 with dipolar moment q, s(r)=− q, ∇δ(r − r 0 ) .
The subtraction method consists in avoiding the singularity in s(r) by separating Ω into two subsets, one surrounding the source, namely Ω 0 , with homogeneous conductivity σ ∞ = σ(r 0 ), and the other being Ω c =Ω\Ω 0 (the complement of Ω 0 in Ω) with electrical conductivity σ c (r)=σ(r) − σ ∞ . This allows to express the electric potential as the sum of two terms, u(r)=u c (r)+u ∞ (r), where u ∞ (r) is the singularity potential generated by a source in an unbounded homogeneous conductor with conductivity σ ∞ (for which analytical expressions exist), and u c (r) is the correction potential satisfying ∇· σ(r)∇u c (r) = −∇ · σ c (r)∇u ∞ (r) (r ∈ Ω), and subject to σ(r)∇u c (r),ň(r) = − σ(r)∇u ∞ (r),ň(r) (r ∈ Γ) [3], [4]. The problem then turns to find the correction potential by means of a numerical method, after which the singularity potential is added.
2) FEM Discretisation: The finite element (FE) formulation of the EEG-FP relies on the variational form of the subtraction version. This can be obtained by multiplying the corresponding differential equation by a test function v belonging to a suitable space H, and then integrating over Ω [3], [4]. After applying the divergence theorem and utilising the boundary condition, the variational formulation results in finding u c (r) ∈ H such that, for all v(r) ∈ H, satisfies a (u c ,v)=l(v), where a : H × H → R is the bilinear form defined as and l : H → R is the linear form given by To proceed with the FEM, a discretisation T of Ω is required. This discretisation is composed by a set of nodes p i (i =1,...,N) and a set of elements t j (j =1,...,N e ) defined upon these nodes. This tessellation allows to construct a discretised FE space V N ⊂ H where to find the numerical solution. More explicitly, we choose V N = span{ϕ i (r): i =1,...,N}, with ϕ i (r) being piecewise functions satisfying ϕ i (p j )=δ ij [3]. Then, we look for u c (r) ∈ V N (an approximation of u c (r) ∈ H) such that, for all v(r) ∈ V N , a( u c ,v)=l(v). This leads to solve the linear system where K ∈ R N ×N is the stiffness matrix with elements K ij = a(ϕ i (r),ϕ j (r)), b ∈ R N is the source vector defined by b i = l(ϕ i (r)), and u c ∈ R N is the vector containing the numerical approximation of the correction potential on the mesh nodes (i.e., u c (r) ≈ u c (r)= N i=1 ϕ i (r)u c i ). The most common scenario available in the literature (and the one adopted in this work) is to utilise linear basis functions ϕ i (r) (i =1,...,N) defined over a tetrahedral discretisation T .Inthis case, the stiffness matrix K can be easily found by assembling element matrices that are computed without the need of numerical integration schemes [12]. However, this is not the case for b, since it depends on the non-linear function ∇u ∞ (r).T os o l v e this issue, Drechsler et al. [4] proposed the use of a Gauss-Jacobi integration scheme, leading to the so called full subtraction (FS) approach. Although they showed great levels of accuracy, the use of a numerical integration algorithm was found to be a timeconsuming process. This makes the FS method unsuitable for real-case scenarios where detailed models composed by tens of millions of elements and houndred of thousands of sources are required [8], [13]. This limitation was tackled by Wolters et al. [3], who proposed to replace u ∞ (r) with π h u ∞ (r),t h e projection of u ∞ (r) in the FE space. This approximation, named projected subtraction (PS) approach, presented similar accuracy than the FS method for high quality meshes, while reducing the computation time notably [4].

3) Sensitivity of the PS Approach to Mesh Quality:
The computation of both K and b in (3) is performed by approximating the integrals (1) and (2) over the corresponding simplices. For clarity, we split the source vector b into two terms, one corresponding to the volume integral in (2), namely b v , and the other to the surface integral, b s . Then, the discretisation is done by computing the element matrices K e ∈ R 4×4 and b v e ∈ R 4 for each tetrahedron t in T , and b s e ∈ R 3 for each surface triangle. The global linear system is finally obtained by properly adding these element arrays [14].
The geometry of the mesh T will have different effects on each of these matrices, and therefore on the numerical solution. In the case of the stiffness matrix, poorly shaped tetrahedra will affect its condition number. This number is an indicator of the problems that we may encounter when solving the global system. Large condition numbers will imply large roundoff errors when using a direct method to solve the global system, or slower performance in case of utilising iterative solvers [1]. The latter are generally preferred in the field for requiring less memory than the former, which may become prohibitively expensive for highly refined models with a large number of degrees of freedom. Although the condition number is one of the limiting factors in the accuracy of the EEG-FP, its effects can be reduced, in the case of indirect methods, by properly preconditioning the system [6]. Moreover, some Krylov subspace (iterative) methods, as the conjugate gradient method, can handle ill conditioned systems provided that this issue is generated by few bad elements [1]. Then, the effect of degenerated tetrahedra strongly depends on the method used to solve the linear system of equations.
The errors associated with the increased condition number of the stiffness matrix will be shared between both PS and FS approaches (as well as any other FEM formulation). However, there is another error that is related to the PS method only. In this approach, we approximate u ∞ (r) with π h u ∞ (r), and consequently ∇u ∞ (r) with ∇π h u ∞ (r). This apparently safe interpolation is strongly linked to the mesh, and impacts on both b v e and b s e , and consequently on b. To show this, we first need to find the relation between the errors in the element vectors and the interpolation error, i.e., the error obtained when approximating ∇u ∞ (r) with ∇u ∞ (r). For clarity, we analyse the impact on both element source vectors separately, starting with the volume integral. The following Lemma (shown in Appendix A) summarises such relation.
denote the relative error of the volume element vector approximation b v e obtained by replacing ∇u ∞ (r) with ∇u ∞ (r). Then, for every t ∈ T , there exists r * ∈ t such that where we defined g(r) ∞ = max r∈t g(r) 2 for any g : is the largest singular value of σ c , and C ϕ is a constant depending on the basis functions within t. Lemma 1 states that the relative error in the volume element vector is upper bounded by the error corresponding to the approximation of the gradient of the singularity potential. This apparently impractical result becomes relevant when relating the errors in such approximation with the geometrical characteristics of T . This was elegantly done by Shewchuk [1], who showed that, in the case of using the PS approach (i.e., ∇u ∞ (r)=∇π h u ∞ (r)), the maximum error within any tetrahedron satisfies where l ij is the length of the edge between nodes i and j (1 ≤ i<j≤ 4), A i is the area of the face opposite to node i (i = 1,...,4), V is the volume of the tetrahedron, and c t is a bound on the curvature constraint of u ∞ (r) in t.
A close analysis to the inequality (5) shows that the upper bound on RE( b v e ) depends on the shape of the element. For a regular (i.e., equilateral) tetrahedron with edge-length a, A = a 2 √ 3/4 and V = a 3 √ 6/12, leading to RE( b v e ) ≤ C t a, with C t constant for each tetrahedron. This suggests that, for perfectly-shaped tetrahedra, this error can be reduced simply by minimising the size of the elements. Unfortunately, this is not always the case. Badly shaped tetrahedra may result in this bound tending to infinity. This can be achieved, for example, in the case of slivers, characterised by V → 0 with no face following that tendency [11].
The impact of the approximation of the singularity potential on b s e can be studied in a similar way. Lemma 2 describes the relation between the relative error in the surface element vector and the interpolation error.
Lemma 2: Let T s be the set of triangles defining the boundary of T . Also, let RE( b s e ) denote the relative error of the surface element vector approximation b s e obtained by replacing ∇u ∞ (r) with ∇u ∞ (r). Then, for every triangle t s ∈ T s , there exists r * ∈ t s such that where C θ (r * ) is a constant depending on the basis functions set over t s . Lemma 2 can be exploited to understand the relation between the RE and the geometry of T s . In case of considering the PS approach, Shewchuk [1] showed that the interpolation error over triangles satisfies where l max , l med , and l min are the maximum, median, and minimum edge lengths of the element, respectively. The bound (7) shows that the relative error on b s e is extremely sensitive to the mesh geometry. In case of generating a surface mesh composed by equilateral triangles with edge length a, A = a 2 √ 3/4, and then RE( b s e ) ≤ C s a, with C s constant for each triangle. However, deformed triangles may impact on the RE differently. This can be achieved, for example, by an isosceles triangle with one angle near π radians, for which A tends to zero while the sides change little.

1) Motivation:
The errors introduced by the PS approach can be mitigated by interpolating the gradient of u ∞ (r) instead of the potential itself. In other words, we compute analytically ∇u ∞ (r) in the mesh nodes, and then approximate its value within the element with the linear interpolant ∇u ∞ (r)=π h ∇u ∞ (r). The reason for this selection is based on the fact that the error bound of a linear interpolation function over a simplex is less sensitive to the element's shape than the corresponding to the gradient of the interpolant [1]. This is valid for both volume and surface element vectors. In the case of b v e , the upper-bound on the error turns out to be where r mc is the radius of the min-containment sphere of the element (i.e., the smallest sphere enclosing the tetrahedron), and c * t is a bound on the directional curvature of the partial derivatives of u ∞ (r) in t (see Appendix B).
The bound given in (8) is far simpler and more robust to the element's shape than the corresponding to the PS. This is evidenced by the absence of the volume, area, or edge-length variables in its denominator. Unlike the bound (5), the error is limited by the element's size only, which can be arbitrarily reduced without the need of specialised mesh-refinement tools. In the particular case of considering an equilateral tetrahedron with edge-length a, r mc = a √ 6/4, and the upper-bound on the interpolation error (and consequently on b v e ) reduces to C * t a 2 , with C * t constant for each element. This represents a clear advantage over the PS method for perfectly-shaped simplices.
In case of interpolating over triangles, the error bound is Again, this is a great improvement over the PS approach, which can be appreciated on the finiteness and independence of the upper bound on any mesh feature other than its size. Moreover, the rate of convergence of RE( b s e ) is increased to a 2 for equilateral triangles, instead of a as found for the PS method.
2) Mathematical Formulation: The discretisation of the PGS approach is obtained by replacing ∇u ∞ (r) with N j =1 ϕ j (r)∇u ∞ (r j ) in the source vector. After some algebraic manipulations, we get where and the vectors g (k ) ∈ R N (k =1, 2, 3) have elements g The matrices L (k ) and R (k ) in (10) do not depend on any source parameter. If sources are assumed to belong to the same homogeneous compartment, L (k ) and R (k ) are constant, and therefore can be computed just once. This is usually the case in the EEG-FP, where sources are assumed to be located in the cortex, which is accepted to be electrically isotropic. Since the vectors g (k ) are computed analytically, the computational load of the PGS method is on the same order to that corresponding to the PS approach [in fact, slightly more than three times according to (10)].

C. Experiments
A set of experiments were performed to study the differences in accuracy and stability between the PS and the PGS approaches. These were separated into local or global experiments, in which we analysed differences in the element and Fig. 1. Scheme of the tetrahedron deformation method used in the second local experiment. An originally-equilateral tetrahedron formed by nodes p 1 , p 2 , p 3 , and p e 4 is deformed by rotating one node (p e 4 ) in θ radians with respect to the circumcentre of the original equilateral tetrahedron (C) and in the opposite direction to one of the fixed nodes (p 2 ). A sliver is obtained if θ takes its maximum value (i.e., θ max = π/2+ sin −1 (1/3), leading to nodes p 1 , p 2 , p 3 , and p s 4 ). The circumradius is the same for any rotation angle, and the radius-edge ratio takes a maximum value of 0.8018 (achieved in the case of a sliver).
global source vectors, respectively. The experiments described below are focussed on exploring the impact that degenerated tetrahedra have on the volume source vector, avoiding further analysis of the surface term. The reason for this is that, unlike volumetric meshes, triangular surfaces can be generated in such a way that degenerated elements are avoided. This can be done, for example, by applying node repulsion algorithms to the mesh [9].

1) Local Effects:
We first tested the influence of the size of the tetrahedra in the accuracy of b v e . This is of great importance for comparing the requirements on element sizes needed by both PS and PGS approaches to achieve a certain relative error. To do so, we considered a source located in the origin with moment [10, 0, 0] nAm, and assumed homogeneous and isotropic electrical conductivities σ c and σ ∞ . Then, we calculated the maximum side-length of an equilateral tetrahedron needed to achieve RE( b v e ) ≤ 0.01, which is independent of the electrical conductivity values. Test elements were placed in such a way that their centroids belonged to the z =0plane, with x and y in the range [0, 0.09] m. This covers a range of distances from the source to the electrode position as found in experimental conditions. The relative errors were then computed for the PGS and PS approaches considering the solution obtained with the FS method as a reference.
Secondly, we studied the influence of the element shape in the computation of the element volume source vector by means of both PGS and PS techniques. Tetrahedra with different shapes were used to compute the relative error with respect to the result obtained with the FS approach. To generate them, we started with an equilateral tetrahedron with side-length a = 0.001 m. Then, we degenerated it by moving one of the vertices towards the plane containing the other three, which were fixed to the initial positions (see Fig. 1). This movement was characterised by a rotation in θ radians with respect to the original circumcentre (C) in the opposite direction to one vertex (p 2 in Fig. 1). This allowed us to generate tetrahedra ranging from equilateral (θ min =0)t os l i v e r s( θ max = π/2+sin −1 (1/3)). The moving node was constrained to the circumsphere of the original equilateral tetrahedron, allowing to generate elements with similar radius-edge ratios, defined as the ratio between the circumradius and the shortest edge length. This is important since this factor is utilised by some mesh-generation algorithms in the field as an indicator of mesh quality (e.g., [10]). In the present case, the circumradius was constant for all elements and equal to a √ 6/4. The minimum side length ranged between a 7/12 (obtained for a sliver) and a √ 6/4 (in case of an equilateral tetrahedron). This led the radius-edge ratio to be constrained in the range √ 6/4, 3 √ 14/14 , which is more than acceptable compared to the ratios reported in the literature (all above 1, see [4], [6]). Elements with different shapes were then translated to positions in the xy plane as done in the previous experiment.
2) Global Effects : We analysed the performance of the subtraction methods in multi-layered spherical head models. Although simple, these representations have the advantage of counting with analytical solutions with which to compare [2]. We modelled the head as a multi-layered sphere representing the scalp, skull, cerebrospinal fluid (CSF), and brain (grey and white matter). The outer radii were 0.092 m, 0.086 m, 0.08 m, and 0.078 m, respectively. The electrical conductivities were considered isotropic for the scalp, CSF, and brain compartments, and set to 0.33 S/m, 1.79 S/m, and 0.33 S/m, respectively. The skull was modelled as an anisotropic layer with radial/tangential conductivities of 0.0042/0.042 S/m. These values were selected from the relevant literature [15]- [18].
Triangular surface meshes were generated using the Distmesh toolbox [9]. They included 144 nodes representing the sensing positions, which were uniformly distributed on the outermost surface using an analytically exact spiral scheme [19]. The surface tessellations were then used to generate the corresponding tetrahedral models using the ISO2Mesh toolbox [20], a Matlab wrapper of Tetgen [10]. Meshes were built to achieve a maximum radius-edge factor of 1.2. We considered a coarser mesh resolution in the brain layer since b v vanishes in it for sharing the same electrical conductivity as the source neighbourhood [4]. Seven models were built, with the number of nodes between 39,000 and 950,000. Each model was generated by properly refining the surface and volume discretisation parameters without any local refinement [4] (Fig. 2).
The mesh quality was assessed using the normalised aspect ratio, defined as q(t)= √ 3h min / √ 2l max , with h min being the shortest height of the tetrahedron [11]. The closer the quality of an element is to the unity, the better shaped the element is. For the case of the deforming tetrahedron in Fig. 1, a simple calculation shows that q(θ) ≈ 1 − θ/θ max . Degenerated tetrahedra were defined as those with q(t) ≤ 0.1. Table I summarises the resulting models.
We used these tessellations to simulate highly eccentric dipolar sources located in the innermost layer. These sources are known to generate the largest errors, and therefore a good indicator of the performance of the technique [3], [4]. We simulated 50 radially-and tangentially-oriented sources uniformly distributed on a sphere with radius 0.0743 m, resulting in an eccentricity of 0.9524. The amplitude was assumed as 10 nAm.
Results were used to compare the subtraction methods as a function of the model refinement. Such comparisons were done in both electric potential and volume source vector. In the first case, we computed the potential for the three approaches and evaluated the relative error utilising the analytical solution as the reference [2]. In the second experiment, we compared the error contribution of different element volume vectors to the global volume source vector as a function of the element quality, utilising Model 5. This was done by calculating the absolute error on the element vectors for both PS and PGS approaches, considering the result obtained with the FS as the reference. Then, we computed the average error introduced by elements belonging to the quality intervals Q i =((i − 1)/10,i/10] (i = 1,...,10). This result was then used to compare the impact of different element qualities to the overall source vector, and consequently to the resulting electric potential.

3) Illustration in a Realistic Scenario:
We applied the method developed to illustrate the impact of the approach used to solve the EEG-FP in the solution of the EEG-IP. To this end, we utilised evoked response potentials from five healthy individuals triggered by the presentation of "happy and angry" faces. The response, known as N170, is a cortical marker specifically linked to facial processing, with neural generators in the fusiform gyrus and superior temporal sulcus [21]. EEG signals were recorded at 500 Hz with a Biosemi 128-channel Active Two system. The experimental setup, data acquisition and preprocessing protocols were based on our previous studies in the field [21], [22]. A detailed head model was built based on the ICBM 2009 atlas [23]. A mesh with approximately 7.9 million tetrahedral elements was created using ISO2Mesh, of which 22,497 resulted degenerated. Isotropic conductivities were assumed for the 7 tissues included in the model: skin (0.435 S/m), fat (0.078 S/m), bone (0.0064 S/m), marrow (0.0286 S/m), cerebrospinal fluid (CSF; 1.79 S/m), grey matter (GM; 0.33 S/m), and white matter (WM; 0.142 S/m). As before, the electric conductivity values were selected from the relevant literature. A slice of the head model is shown in Fig. 7(a).
Lead field matrices were computed utilising both PS and PGS methods for 33,255 dipolar sources located on a surface in between the GM/CSF and WM/GM interfaces, and unconstrained in orientation [24]. The source strength was 10 nAm. The EEG-IP was solved utilising the standardised, low-resolution, brain electromagnetic tomography algorithm (sLORETA) [25] considering both lead field matrices.

4) Implementation:
We implemented the FE formulations in Matlab 2015a (The MathWorks, Inc., Natick, Massachusetts, United States). The computation of the numerical integrals involved in the FS approach was based on the Gauss-Jacobi method. To this end, we used the jacpts function form Chebfun [26] for getting the Gauss-Jacobi quadrature nodes and weights of arbitrary order (we considered an integration order equal to 4 in the experiments). Linear systems were solved using the preconditioned conjugate gradients method with a tolerance of 10 −10 and with incomplete LU preconditioning. All experiments converged to the result.

A. Local Effects
The maximum edge length needed for an equilateral tetrahedron to achieve a relative error less or equal to 0.01 is shown in Fig. 3. It can be appreciated that the PS method needs proper mesh refinement for reaching such error bound. This was found the case for almost any test point, apart from those in which ∇u ∞ (r) is practically linear, and therefore, well approximated by ∇π h u ∞ (r) [yellow bands in Fig. 3(a)]. On the other hand, the PGS approach is almost insensitive to the element size for tetrahedra located at a distance greater than 3.7 cm to the source. This means that the PS technique will require smaller elements than the PGS approach for achieving a certain RE, highlighting the increased accuracy provided by the latter. Consequently, the PGS method allows the user to reduce the element density in regions not belonging to the source neighbourhood without compromising the results. This is not the case for the PS approach, which will clearly benefit from a fully-refined model. Less than 1% of the test points were found to have a lower edge length in the case of utilising the PS method. Fig. 4 shows the resulting maximum relative error on the computation of b v e for a given distance to the source as a function of the shape of the element, for both PS and PGS methods. It can be seen that, as expected, the element source vector computed with the PS approach is unstable for degenerated elements (i.e., θ → θ max ). On the other hand, the PGS technique is not only insensitive to the element's shape, but also more accurate than the PS method in at least 2 orders of magnitude. Fig. 5 presents the relative errors of the numerical solutions of the EEG-FP obtained using the three subtraction approaches. Results are presented for the seven models described in Table I. As expected, errors decrease at different rates for the three methods, being the FS approach the one leading to more accurate results. On the other hand, the PS method provides solutions with RE greater than 10% for every model as a consequence of the presence of non-regular tetrahedra in the mesh. This rate is hugely outperformed by the PGS technique, which converges to the results obtained with the FS method despite of the number of badly-shaped elements. Regarding the computa-tional effort needed by each method, the FS approach required (in average for the whole experiment) 180.72 times the effort needed for the PS method, whereas the PGS approach only required 3.23. Similar conclusions can be extracted from the normalised relative difference measure (NRDM) and magnification (MAG) errors, presented in the Supplementary document.

B. Global Effects
The histogram representing the mean absolute error on the element volume vector as a function of the mesh quality interval is presented in Fig. 6. It can be noted that, in the case of the PS method, the error is monotonically decreasing, indicating that the mesh quality has a huge impact on the assembled volume source vector. This is not the case of the PGS approach, which exhibits a nearly flat error for elements belonging to the interval (Q 1 ,Q 6 ).

C. Realistic Scenario
Solutions to the EEG-IP utilising the lead field matrices computed with the PS and PGS methods are presented in Fig. 7(b) and (c), respectively. For both approaches, sources were found mostly in the right fusiform area, as previously reported [22]. However, differences in the estimated standardised current density power maps can be appreciated, mostly related to MAG errors in the PS technique. Although such differences were not found significant in the source localisation results, they resulted important for properly estimating the source strength.

IV. DISCUSSION
The analysis presented here demonstrates the instability of the PS approach with respect to the mesh geometry for computing the source vector in the EEG-FP. More specifically, we derived upper bounds on the RE of the element source vectors, which were found to tend to infinity for highly degenerated elements, such as slivers. The reasons are based on the unstable nature of linear interpolants over deformed simplices, as thoroughly investigated by Shewchuk [1]. The value of the upper bounds for characterising interpolation errors depends exclusively on their tightness; the tighter the bound is, the better it will represent such error. In the present case, upper bounds (7) and (9) are known to be tight to within a factor of three, and bounds (5) and (8) were suggested to share such characteristic [1]. This allows us to conjecture that the bounds introduced by Lemmas 1  and 2 are equally tight. This hypothesis was confirmed by computer simulations, which showed that the RE on the element volume vector was highly influenced by the element's shape, as predicted by the Lemmas. The impact of degenerated elements was then found to have a detrimental effect on the global source vector (see Fig. 6), and consequently in the solution of the EEG-FP (see Fig. 5).
To mitigate these issues, we presented the PGS approach, in which we proposed to project the gradient of u ∞ (r) onto the FE space. This was based on the finiteness of the corresponding error bounds for this case, which are independent of the element's geometry [see (8) and (9)]. A comparative analysis between the PS and PGS methodologies for computing the element volume vectors showed that the latter was not only less sensitive to the shape of the tetrahedra (see Figs. 4 and 6), but also more accurate in the case of assuming regular elements (see Fig. 3). This is grounded on the fact that the PGS technique performs a linear interpolation of ∇u ∞ (r) over each simplex, whereas the PS technique assumes it constant. Such difference is also evidenced in the upper bounds obtained for regular elements, converging with rate a 2 in case of using the PGS method, rather than a as obtained for the PS approach.
The benefits of the PGS methodology resulted in an increased accuracy on the solution of the EEG-FP with respect to the PS approach, as presented in Fig. 5. These benefits were more evident as the number of nodes became larger, independently of the quality of the mesh. Such characteristic was confirmed in Fig. 6, were we found almost no difference in the mean error introduced by the PGS technique for elements in the range (Q 1 ,Q 6 ).O n the contrary, the PS approach showed a pronounced decrease of this error as a function of the element quality, indicating its sensitivity to low quality tetrahedra. Similar results were obtained for all the models described in Table I, indicating the robustness of the proposed technique to the mesh quality.
We demonstrated that the errors introduced by low quality elements are detrimental for the electric potential solution computed with the PS approach. Although degenerated simplices with q(t) ≤ 0.1 introduced the highest errors, low quality elements (but not necessarily degenerated) were shown to impact on the solution as well. This is of utterly significance since  low quality tetrahedra are produced by the most popular mesh generators, as those based on the Delaunay triangulation algorithm, either in their classical and constrained versions [9]- [11]. Although there exist several algorithms and methodologies for generating and/or improving the quality of tetrahedral meshes (e.g., Cleaver [27], Stellar [28], Distmesh [9], NETGEN [29]), it is not yet available a technique that corrects all degeneracies in a robust and consistent manner [11]. This means that results will vary depending on the utilised method, which can even introduce further degeneracies in the process. Nevertheless, these algorithms were successfully used for obtaining highquality meshes with which to solve the EEG-FP by means of the PS approach [4], [6]. However, it is not a standard practice in regular EEG applications to apply mesh quality boosting algorithms before running the corresponding experiments. The reasons are that such methodologies are not straightforward and, more importantly, not properly acknowledged (or even utilised) in the literature. In these regards, the PGS method presents a viable, efficient, and robust alternative to compute lead field matrices without increasing the computational complexity.
Results from the idealised experiments utilising spherical head models (see Fig. 5) suggest that the use of the PGS approach in a mesh discretisation composed by 440,000 or more nodes should be enough for reaching an almost negligible error in the resulting lead field matrix, and consequently in the EEG-IP (less than 5 mm localisation error according to [16,Fig. 6]). However, multiple factors other than the algorithm used for solving the EEG-FP are known to impact on the calculation of the lead field matrix in realistic scenarios, each of them in a particular way. These factors include (but are not limited to) the quality and size of the discretisation (both minimised by the PGS method), the consideration of a fully-realistic head and electrode model [13], [30], the uncertainty in the electrical conductivity field [16], [17] and electrodes' position [31], and the source model and location [32]. These errors will consequently affect the solution of the EEG-IP, in addition to those introduced by the pre-processing workflow and inverse algorithm used. Therefore, efforts are needed to shed light on the relative importance of different factors in source localisation, as well as in understanding the complex interplay between both forward and inverse problems.
It is very important to point out that the improvements provided by the PGS method come at the expense of very little extra computational effort. Such characteristic is essential in realistic scenarios, in which over a million source vectors need to be calculated. This is not the case of the FS approach, whose computational requirements may become prohibitive. We showed that the PGS technique delivers fast and accurate results, making it a truly competitive method in the field. We believe that the approach presented here will help in repositioning the subtraction method as an efficient technique able to produce high quality lead field matrices even in the presence of degenerated elements. Further work will include the analysis of the impact of the mesh quality in the computation of the return currents, which are a secondary source of signal in magnetoencephalography.

V. C ONCLUSION
We have shown that numerical solutions of the EEG-FP computed with the PS method are highly sensitive to low quality tetrahedra. The reasons are based on the instability of the error bounds on the element source vector when the gradients of linear interpolants are computed. To solve this problem, we presented the PGS approach, in which we project and interpolate the gradient of the singularity potential onto the FE space. This selection was based on the stability of the corresponding error bounds, which were found to be independent of the mesh geometry, and therefore stable regardless of the mesh quality. Analytical results and in silico experiments allowed us to show the advantages of the PGS approach over the PS method, which resulted in a better performance even under high-quality tessellations.

APPENDIX A DERIVATION OF LEMMA 1
To derive an upper bound of RE( b v e ), we work separately with its numerator and denominator. In the case of the numerator, we have where θ i (r) is the angle between both vectors within the inner product. Applying the inequality Ac 2 ≤ A 2 c 2 = λ max (A) c 2 in (13), valid for any A ∈ R N ×N and c ∈ R N [33], and noting that g(r) 2 ≤ g(r) ∞ for any vector function g : Finally, using the Cauchy-Schwarz inequality, we obtain Ω cos(θ i (r))dr 2 ≤ V Ω cos 2 (θ i (r))dr ≤ V 2 , and then b In the case of the denominator, the integral version of the mean value theorem allows us to say that there exist r * ∈ Ω e such that Ω e f (r)dr = Vf(r * ). Then, It is clearly seen that, since the cosine terms cannot be simultaneously zero, C ϕ is upper bounded. Lemma 2 is shown similarly, with C θ = | cos(α(r * ))| −1 , and α being the angle between σ ∞ ∇u ∞ (r * ) andň(r * ).

APPENDIX B NOTES ON THE BOUNDS
In case of approximating ∇u ∞ (r) with ∇u ∞ (r)= ∇π h u ∞ (r) (as in the PS approach), expressions on the right hand-side column of [1, Table 2] can be utilised without mod-