One Definition to Join Them All: The N-Spherical Solution for the EEG Lead Field

Albeit its simplicity, the concentric spheres head model is widely used in EEG. The reason behind this is its simple mathematical definition, which allows for the calculation of lead fields with negligible computational cost, for example, for iterative approaches. Nevertheless, the literature shows contradictory formulations for the electrical solution of this head model. In this work, we study several different definitions for the electrical lead field of a four concentric spheres conduction model, finding that their results are contradictory. A thorough exploration of the mathematics used to build these formulations, provided in the original works, allowed for the identification of errors in some of the formulae, which proved to be the reason for the discrepancies. Moreover, this mathematical review revealed the iterative nature of some of these formulations, which allowed us to develop a formulation to solve the lead field in a head model built from an arbitrary number of concentric, homogeneous, and isotropic spheres.


Introduction
The calculation of the distribution of electric potential due to the current in a conductor is not trivial. However, this calculation is necessary to reconstruct, for example, the cortical sources of activity of electroencephalography (EEG). Currently, the solution of this forward problem is usually calculated using numerical approaches, where the conductor (in this case, the head) is divided into a continuous series of elements of piecewise constant conductivity. However, for very simple cases where the geometry can be mathematically parametrized, it is possible to use a mathematical approach, which would give an "exact" solution for the problem. Since the head can be approximately modeled as a sphere, the widespread solution is that of a spherical conductor. While it is not precise for a realistic head, it is exact for spherical geometry.
The electrical lead field in the surface of a homogeneous sphere was mathematically solved by Fran in 1952 [1], and the solution for a multilayer concentric set of spheres was proposed by Arthur and Geselowitz in 1970 [2] and simplified by Cuffin and Cohen in 1979 [3]. The definition of the geometry-defined constants in the solution was not precise, and it was later improved by Lütkenhöner [4] and included in the FieldTrip toolbox [5]. In parallel, Yao and colleagues [6,7] and Nunez and Srinivasan [8] developed their own approaches. However, these three approaches seem, at first, incompatible. In addition, Lütkenhöner's work is not publicly available, and Nunez and Srinivasan's solution contained several errors, later corrected by Naess [9]. Right now, it is not clear if all these solutions define the same phenomena or if some of them are erroneous.
Although the increase in computational power allows for complex head models for solving the EEG forward problem, the use of computationally simple spherical solutions continues to be broadly extended [10,11]. For example, the spherical solution is still prevalent in cases where the accuracy of the forward solution is not essential [12] or as part of a two-step iterative procedure, where the simplified geometry is used in a first approach and a more realistic one is used in the last steps. Moreover, the existence of accurate analytical solutions is of paramount importance for the evaluation of the accuracy of numerical approaches. As these solutions are naïve to the parametrization of the geometry, their error will be similar in simple and realistic cases, and the concentric spheres model has been broadly used as testing ground [13][14][15][16].
In this work, we go over the three most common formulae for solving the lead field in a concentric spherical head model [3,6,8] and compare the different solutions. In addition, we create an alternative formulation that explodes the iterative nature of the previous solutions and allows for an arbitrary number of concentric spheres. Then, we evaluate the validity of this approach, comparing the results with implementations based in the previous literature. Finally, we provide an implementation of the algorithm to be freely used by the reader.

Mathematical Development
For a point-like current (i.e., an elemental dipole) P, in an infinite medium of electrical conductivity σ, located in the Z-axis and oriented radially, the potential in any point of space → r is given by the solution to Laplace's equation, which in spherical coordinates can be written as [2] where r 0 is the distance from the dipole to the origin, r and θ are the radial and azimuthal coordinates of point → r , and P l is the Legendre polynomial of degree l. As the orientation of the dipole is radial (i.e., with the Z-axis), the potential distribution is rotationally symmetric, and the result is independent of ϕ.
When the medium is replaced with a set of s concentric spheres, of piecewise constant conductivity σ s and surrounded by vacuum (i.e., no conductivity), the solution for the potential is [3] where s denotes the region, and A s l and B s l depend on the boundary conditions, particularly the position of the dipole, the radii of the spheres, and the conductivity of the compartments. As these constants specifically affect the terms of the spherical harmonics' expansion given by P l , they are sometimes referred as a spherical harmonics spectrum or spatial filters [6]. We can relate them with the boundaries between compartments: when the potential generated with a dipole P reaches a boundary between media, the behavior can be described as a double layer of charges distributed in the boundary, with one in the i-th medium and another in the i + 1-th medium [17]. Then, this distribution of charges will themselves generate a perturbation in the electric field that will propagate both inwards and outwards [8], and A s l and B s l describe this propagation. It is important to note that, in the following equations-and in the same fashion as in the definitions of A s l and B s l above-the subscript will represent the degree of the Legendre polynomial (l, in general), and the superscript will represent the boundary where the constant is calculated (s, in general).
Applying the boundary conditions (i.e., continuity of the current and potential across boundaries and no current leaving the outer sphere of the conductor), Cuffin and Cohen [3] provided a closed formula to calculate the potential on the external surface of a four-layer conductor. This formula proved later to be inexact and was corrected by Lütkenhöner [4]. Using a similar approach, Nunez and Srinivasan [8] (later corrected by Naess and col-leagues [9]) developed a formula for the potential at any point in the conductor. The general formulation is similar to (2): In their work, the authors provide formulae to calculate the terms A s l and B s l for a four-layer conductor (see the Appendix in [9] for details on the derivation of (4)- (14)): Here, R s is the radius of sphere s, and σ s is the conductivity of the compartment bounded by it. The constants . Z l , Y l , and V l , for their part, are proportionality constants relating the terms A s l and B s l for s greater than 1.
One of the main caveats of this formulation is that the position of the dipole is included in the terms. Thus, it must be calculated for each dipole position. However, we can rewrite the formulae in (4)-(11) extracting the dipole position as follows: etc. The same idea can be applied to B s l to obtain the terms ∼ B s l . Now, we rewrite the potential in (3) as follows: We can further simplify the definitions for the constants by introducing a new constant, U l , into (12)- (14): This formulation reveals the iterative nature of these constants. For the external compartment, the radius (R 5 ) is infinite, and the conductivity (σ 5 ) is zero. Therefore, the first constant, U l , has the most straightforward definition. The rest of the constants depend on the previous one. However, the innermost constant, Z l , does not follow this rule. We can also observe that, in (16), the definition of .
We observe that the new constant In this case, the innermost constant, ∼ A 1 l , is the simplest, and the other constants depend on the previous one.
In this form, if we take r 0 r l+1 as a common factor, the constants can be seen to match the terms for the spatial filter W s l in [7] ∼ A s l r 2l+1 Finally, taking advantage of the iterative nature of the definitions as mentioned above, we can write the constants for a conductor consisting of n concentric spheres as

Evaluation of the Accuracy and Performance of the Proposed Implementation
In Section 2 of this work, we have developed an iterative formulation for the calculation of the potential at any point of a spherical conductor made from concentric spheres of homogeneous materials. In this section, we evaluate the accuracy of our formulation, compared with the classical formulations by Lutkenhoner, Yao, and Naess, together with the computational cost associated with each formulation. In this evaluation, we compare the different results using Lutkenhöner's constant, as implemented in the FieldTrip toolbox [5], as the gold standard. The code for these comparisons is available in the GitHub repository of the work (https://github.com/rbruna/eeg-spheres/, accessed on 1 June 2023), where we provide implementations for the lambda constant using the mathematics described in the different works referenced in the main body of the manuscript. Figure 1 shows the error (compared to Lutkenhöner) as relative error (difference of values divided by the magnitude of the value) when using the classical values for the conductivities of the brain (1/3 S/m), the skull (1/(3 × 80) S/m), and the scalp (1/3 S/m). The error of all methods is similar and is inside the expected error for double-precision float point arithmetic. the different results using Lutkenhöner's constant, as implemented in the FieldTrip toolbox [5], as the gold standard. The code for these comparisons is available in the GitHub repository of the work (https://github.com/rbruna/eeg-spheres/, accessed on 1 June 2023), where we provide implementations for the lambda constant using the mathematics described in the different works referenced in the main body of the manuscript. Figure 1 shows the error (compared to Lutkenhöner) as relative error (difference of values divided by the magnitude of the value) when using the classical values for the conductivities of the brain (1/3 S/m), the skull (1/(3×80) S/m), and the scalp (1/3 S/m). The error of all methods is similar and is inside the expected error for double-precision float point arithmetic.  [6], Yao (2001) [12], Yao (2003) [7], and Naess (2017) [9], Bruna (2023): this research.
For a more general view, Figure 2 shows the relative error of each definition of the constant when using random (yet symmetrical) values for the conductivity. We can observe that the definitions of Yao in their 2000 and 2001 paper have the lowest relative error, while their definition (corrected as reported in the main body of the manuscript) in the 2003 paper had the largest one. The relative error in our proposed method is similar to that of Naess, which makes sense, as their definition is our main inspiration. Nevertheless, in all the cases, the relative error is inside the numerical error of the double-precision floating point operations.  [6], Yao (2001) [12], Yao (2003) [7], and Naess (2017) [9], Bruna (2023): this research.
For a more general view, Figure 2 shows the relative error of each definition of the constant when using random (yet symmetrical) values for the conductivity. We can observe that the definitions of Yao in their 2000 and 2001 paper have the lowest relative error, while their definition (corrected as reported in the main body of the manuscript) in the 2003 paper had the largest one. The relative error in our proposed method is similar to that of Naess, which makes sense, as their definition is our main inspiration. Nevertheless, in all the cases, the relative error is inside the numerical error of the double-precision floating point operations.  Table 1 shows the times (averaged over 1 million executions) for the calculation of the constant for a three-sphere model. The implementation for Lutkenhöner's definition is a closed formula and is by far the fastest but has the caveat that it can only be used to estimate the potential in the outer interface. The definitions by Yao are next, with the one from 2003 being the fastest (although the one with the larger error, see above). The implementation proposed here is approximately 30% faster than Naess, likely because their closed formula is defined for models with four spheres. In any case, the constant must be generated only once per head model, and the computation time for spherical harmonics of order 100 is less than a millisecond in all cases. Table 1. Average time for the execution, over 1 million runs, of each definition for a three-sphere model with order of the spherical harmonics up to 100.

Discussion
In this work, we propose a reformulation of the solution for the potential in a conductor formed using a set of concentric spheres with piecewise constant conductivity. This reformulation reconciles the different definitions provided by Cuffin and Cohen [3], Nunez and Srinivasan [8], and Yao [7]. In doing so, we have also uncovered the iterative nature of the solution. Therefore, the formulation can be naturally extended to an arbitrary number of spheres.
The original algorithms used as the basis for this work are valid only for a fixed number of spheres, usually three [6] or four [3,7,8]. This limitation makes it necessary to use different algorithms for different models or to trick the algorithm by creating contiguous regions with the same conductivity, wasting precious CPU time. With the proposed for-  [4] and Yao (2000) [6], Yao (2001) [12], Yao (2003) [7], and Naess (2017) [9], Bruna (2023): this research. Table 1 shows the times (averaged over 1 million executions) for the calculation of the constant for a three-sphere model. The implementation for Lutkenhöner's definition is a closed formula and is by far the fastest but has the caveat that it can only be used to estimate the potential in the outer interface. The definitions by Yao are next, with the one from 2003 being the fastest (although the one with the larger error, see above). The implementation proposed here is approximately 30% faster than Naess, likely because their closed formula is defined for models with four spheres. In any case, the constant must be generated only once per head model, and the computation time for spherical harmonics of order 100 is less than a millisecond in all cases. Table 1. Average time for the execution, over 1 million runs, of each definition for a three-sphere model with order of the spherical harmonics up to 100.

Discussion
In this work, we propose a reformulation of the solution for the potential in a conductor formed using a set of concentric spheres with piecewise constant conductivity. This reformulation reconciles the different definitions provided by Cuffin and Cohen [3], Nunez and Srinivasan [8], and Yao [7]. In doing so, we have also uncovered the iterative nature of the solution. Therefore, the formulation can be naturally extended to an arbitrary number of spheres.
The original algorithms used as the basis for this work are valid only for a fixed number of spheres, usually three [6] or four [3,7,8]. This limitation makes it necessary to use different algorithms for different models or to trick the algorithm by creating contiguous regions with the same conductivity, wasting precious CPU time. With the proposed formulation, the algorithm can be naturally extended to the required number of spheres without the need for optimization. In addition, we evaluated the performance and accuracy of the proposed formulation, showing that, while the solution is slightly slower that closed-form formulations, the accuracy is similar to those and close to the accuracy error for floating point arithmetic.
While the iterative nature of the solution allows for the calculation of the potential in a conductor with an arbitrary number of spheres, the practical applications of this feature are negligible. This is because, as the spherical volume conductor for EEG is an acceptable head model for fast calculation of the lead field in a real head, the geometry is not fit to the real head shape; therefore, an increase in the precision of the model will not tend toward an increased precision of the real volume conductor. However, this feature can be exploited in the evaluation of the accuracy of other, e.g., numerical, methods. For example, it is well known that classical BEM models are inaccurate in interfaces between media of different conductivity [18]. Some approaches have been developed to overcome this issue (as the isolated skull approach [19] or the symmetrical BEM methods [16]), but it is not clear how they will behave in more complex scenarios, for example, in the presence of low-conductivity pieces (i.e., fragments of skull). In this vein, the formulation proposed in this work could be of great use to evaluate more complex scenarios, with several high and low conductivities of different thickness layers interleaved.
Finally, in Appendix A, we provide an implementation of the algorithm, programmed in the MATLAB language and compatible with GNU Octave, to be used under the terms of GNU's GPL version 3. This implementation can be compared with the classical formulations, and the code for such a comparison can be obtained from GitHub (https://github.com/rbruna/eeg-spheres/, accessed on 1 June 2023). Appendix B shows the results of these comparisons, both in relative error (compared to Lutkenhöner [4]) and in computation time.