I Introduction
The inverse problem of source reconstruction in electroencephalography (EEG) aims at finding the source distribution that best explains the electric potentials measured noninvasively from electrodes on the scalp surface. A model of the head’s electromagnetic transmission plays a central role in accurate source localization. This model is built from a specification of both the geometry and the conductivity distribution of the modeled tissue compartments (scalp, skull, cerebrospinal fluid, brain gray and white matter, etc.). Because it is impractical to directly measure the head tissues conductivities in vivo for a specific subject, default values are often used. One problem is that the braintoskull conductivity ratio reported in the literature varies from 4 to 80 [1, 2, 3]. As the skull conductivity greatly influences the solution of the forward problem [4], localizing brain sources using uncertain conductivity values leads to important errors [5, 6, 7, 8]. Taking into account the composite structure of the human skull could improve the accuracy of EEG source analysis [9]. In this case, the conductivity of each skull tissue should be estimated simultaneously. A possible solution is to estimate the tissue conductivities and cortical activity simultaneously [10, 11, 12, 13]. Doing this requires to solve EEG forward problem for possibly many different conductivity configurations. A lead field is a linear operator relating the brain electrical activity to the potentials on EEG electrodes. Computing one lead field requires a matrix inversion which is computationally intensive for realistic head geometry represented with meshes with a large number of vertices. Thus, the required time for computing a large number of solutions for different conductivities quickly becomes impractical.
One way to deal with this problem is to approximate lead fields using a relatively small set of precomputed solutions. The reduced basis method approximates the solution of parametrized PDEs [14] by reducing the number of basis functions for Galerkin projection based on a set of already computed exact solutions. Our method is inspired by the reduced basis method, but adapted to the particular structure of EEG forward problem and the nature of its parameter space. In an another method presented in [10]
, only the skull conductivity varies and each element of the lead field matrix is approximated using polynomial interpolation based on a set of precomputed values. With this approach, the complexity of the approximation would increase fast with the number of unknown conductivities.
In this work, we propose a fast lead field approximation method which is robust to the head model complexity and to the number of unknown conductivities.
Ii Forward problem
The useful frequency range for electrophysiological signals in MEG and EEG is typically below 1 kHz, and most studies deal with frequencies between 0.1 and 100 Hz. Consequently, the physics of MEG and EEG can be described by the quasistatic approximation of Maxwell’s equations [15]. In a conductive environment, it yields the fundamental Poisson equation with boundary condition :
(1) 
where is the head domain,
its boundary with outward pointing normal vector
, is the conductivity, is the primary current source density, supposed to be known in the forward problem, is the unknown electric potential. It is usually assumed that the head is composed of several subdomains with different electrical conductivities . We assume a piecewise constant head model (Fig. 1). Conductivity is supposed isotropic and constant in each domain . The different tissue conductivities are denoted .Iia Numerical approximation of the forward problem.
Equation (1) has no analytic solution for realistic (nonspherical) head geometries but several numerical methods can be used to approximate its solution. We will focus on the finite element method (FEM) [17, 18] and the boundary element method (BEM)[19, 16], which are both based on Galerkin projection. The method presented in this paper applies both to FEM and BEM but will be illustrated on the symmetric BEM [16].
For both the FEM and BEM approaches, the initial continuous problem (1) is discretized to a linear system of type:
(2) 
Let be the number of unknowns in the head model. Vector represents unknowns of the model, which are values of the potential on the mesh nodes for FEM and standard BEM, and potentials and their normal derivatives on the meshes for symmetric BEM. The matrix , called ”head matrix”, can be computed once the head geometry, its conductivity values and finite element basis functions are fixed. Vector depends on source configuration. Notice that depends on conductivity. But an important property of matrix is that it can be represented as a linear combination of conductivity independent matrices:
(3) 
where are scalar functions, represents the number of such conductivityindependent components . For FEM . For symmetric BEM multipliers have more complex structure, for example .
IiB Source model
The most common source model to represent electrical activity in the brain is a ”current dipole”. It represents an oriented source of current located at a single position
, with dipolar moment
, and it is denoted by , where is a Dirac delta distribution.The source space can be seen as a finite set of dipoles with known positions. Without loss of generality, we can also assume that the dipole orientation is known. Let represent the amplitude of th source and the th source with unit amplitude:
The source term in (2) is linear with respect to , therefore:
where th column of matrix corresponds to the th unit source and is the vector of source amplitudes. Let us notice that in the case of FEM, matrix does not depend on conductivities. For BEM, it does, so we will note it for generality. Moreover, can be represented as follows:
(4) 
where matrices are independent of and are scalars. In the case of symmetric BEM are .
The linear system to solve then becomes:
(5) 
IiC Lead field matrix
In the context of EEG, the lead field is a linear operator that maps source activation to potentials at sensor locations:
Computing requires applying to a matrix which selects or interpolates potentials only at electrode positions: . Let be the number of EEG electrodes.
Using (5) and the selection matrix :
Let us remark that since the electric potential is only defined up to a constant, head matrix is not full rank and has a onedimensional kernel. So the inverse notation actually implies a deflation which is usually applied to this type of situation.
The lead field matrix can thus be computed as:
(6) 
Each column of the lead field matrix represents the contribution of the corresponding unit norm source on the EEG electrodes.
Iii Fast lead field approximation method
Once and have been computed, the numerical complexity of lead field computation is essentially due to the inversion of the head matrix (see (6)). For realistic head models implying a large number of unknowns the time required for computing lead fields for a large set of conductivity configurations becomes impractical. In this context, it can be worthwhile to compute not the exact lead fields for all needed conductivity configurations, but an approximation thereof.
The general idea of our method is the following. First, we choose a domain of interest in conductivity space, i.e. a range of values for the head tissue conductivities. We then select basis points within this domain, for which exact forward solutions will be computed (Fig. 2). Based on these exact solutions a lead field approximation will be computed for any other conductivity configuration.
In this section, we develop an algorithm which:

selects basis points within the domain of interest;

based on these basis points, computes a fast approximation of the lead field matrix with controlled error.
Iiia Approximation problem formulation
Assume that the approximation of the lead field matrix is parametrized with a vector where is the number of basis points: We denote it . As the lead field is a linear operator, we consider the following approximation error, using the Frobenius matrix norm :
(7) 
This error involves matrix which is not known a priori. We therefore introduce the definition of an upper bound for the approximation error. The idea is to bound error , that we can not directly compute, with an error function , which can be easily computed, converges to zero with increasing number of basis points and thus ensures the convergence of to zero.
Definition 1
For , function is an upper bound approximation of if:
Let us remark that is a constant which does not depend on but may depend on .
Taking into account that is positive, it can be directly derived from definition 1 that if then .
Definition 2
Matrix is a valid approximation of with optimal approximation parameter if both following conditions hold:

[(i)]

,

.
IiiB Choice of the upper bound approximation
As mentioned previously, the point of the approximation is to avoid inversion of the matrix . One possible way to parametrize the approximation matrix is to represent the inverse of the head matrix as a linear combination of precomputed inversions at basis points:
(8) 
Based on (6), the proposed lead field approximation takes the form:
Using the linear decomposition of (4), we represent the lead field approximation as a linear combination of matrices which are independent of conductivity:
(9) 
where and .
As are known, the question is how to compute optimal parameters . Based on the property of an inverse matrix: and on (8), we could seek by minimizing the following expression:
where
is the identity matrix. But in the context of the EEG forward problem, it is more relevant to consider matrix
, which has smaller dimension than and which contains only the relevant part of . Let us introduce following optimization problem:(10) 
The function is a candidate to be an upper bound approximation. Proposition 1 shows that it is indeed the case under a certain condition.
Proposition 1
If matrices are linearly independent, then
is an upper bound approximation of error (7) and is an optimal approximation parameter.
Proof:
Now we can see that is an upper bound of error (7):
Let represent the vectorization of matrix . Using the fact that , problem (10) can be reformulated as:
which is a linear projection problem. It means that if
is a set of linearly independent vectors for
, then . This directly implies that so both conditions of Definition 2 are verified.A linear dependence of matrices would mean that lead fields corresponding to some basis points could be represented as a linear combination of other precomputed lead fields. So the linear independence condition of Proposition 1 depends on selection of basis points. This will be discussed in detail in the next section.
We now express the problem (10) as a simple least squares problem. Using the linear decomposition of as a sum of conductivityindependent matrices (3), we get:
We can then vectorize matrices and ”store” them as columns of a new matrix . Let denote the th column of . The matrix will thus have columns and rows.
Let us also denote , where and is an identity matrix of size :
The following proposition allows to compute :
Proposition 2
The solution of problem (10) is given by the solution of the system:
(11) 
Proof:
Using the decomposition of matrix (3), we get:
Remarks:

The size of matrix is . The complexity of problem (11) therefore does not depend on the size of the head model but only on the number of basis points.

Matrices and are independent of so they can be precomputed. Moreover, adding new basis point amounts to adding new elements to these matrices  they do not have to be fully recomputed.
IiiC Choice of basis points
The method introduced in section IIIB supposes to dispose of a set of basis points with precomputed matrices and . This section answers the question of how to select them. To work with not continuous but discretized domain of interest, we start by sampling it. We denote the set of basis points in conductivity space, which we denote . The idea is to start with a small number of basis points and to add new points in an iterative manner. We propose Algorithm 1.
Algorithm 1 starts by sampling the domain of interest and choosing initial basis points from these samples. One possible choice is to use a dense uniform grid for sampling and to take boundary vertices for initial basis points. Then, for each sample, we compute a valid upper bound error . We choose the sample which results in highest upper bound error as a new basis point and update matrices and for this conductivity configuration. This means that the new basis point is the conductivity sample with the highest error when approximated with the current basis. We iterate this procedure until we reach some convergence criterion. Numerical results show (see Section IVB) that the valid upper bound error decreases with the same speed as the approximation error when increasing the number of basis points. Controlling the upper bound error is a reasonable empirical choice for convergence criterion.
This algorithm has interesting properties which make it more efficient than just a random selection of basis points.
Algorithm 1 controls, in a strong way, the valid upper bound error:
This is true because the projection space at step is a subspace of the one at step .
To ensure that the system 11 has a unique solution we need the matrix to be full rank. Let us analyze two possible situations:

Let Algorithm 1 verify . It means that each iteration adds a basis point which has strictly positive projection error, which means that it is linearly independent of current basis. This implies that the updated basis stays full rank.

If such that , it means that lead fields for every sample from conductivity space can be exactly represented by linear combination of basis elements. In practice, this naturally means that we don’t need to continue adding basis points, unless we introduce new conductivity samples.
All these observations show that Algorithm 1 guarantees that we will verify necessary condition for Proposition 1 and that the system 11 has a unique solution.
Remark: The fact that projection matrix is always full rank, ensures that matrix is always invertible, thus is continuous with respect to conductivity. So is the upper bound approximation . This guarantees that this error is bounded on any compact conductivity domain. And we still need the condition to ensure the existence of the matrix for symmetric BEM.
Once the basis points are chosen and intermediate matrices are computed, a new lead field matrix corresponding to any conductivity in the domain can be approximated with Algorithm 2.
Iv Numerical results
Iva Data description
For numerical implementation of the proposed method, we use the data from [20], which includes anatomical data as well as EEG data recorded with Yokogawa/KIT, and processed with Brainstorm [21]. We use a realistic head model with three layers  brain, skull and scalp. Each surface is represented as a triangular mesh with 1082 vertices. The source space contains 15002 dipoles with fixed orientations (normal to surface).
We use the symmetric BEM implementation from OpenMEEG [16, 22]. The size of head matrix is . Matrix has size .
We use a model with 41 EEG electrodes, so the size of the selection matrix is and the lead field is a matrix.
The conductivity of the scalp is taken equal to 1, while the brain and skull conductivities ( and respectively) are variable and form the parameter space (a subspace of ). We are interested in a subset of this space represented by a rectangle which covers values . This rectangle is uniformly sampled with points.
IvB Approximation error convergence
For each of the 625 sample points, the exact lead field is computed in order to evaluate the approximation error (7). This is just done for purposes of validation and not necessary for the algorithm. We then compare it with its upper bound error, computed with (10).
We follow Algorithm 1 initializing the basis point set with the four corners of the parameter space, and new points are added one by one.
As shown in Fig. (b)b, the approximation error (7) over the conductivity sample points is bounded by with 22 basis points. Fig. (a)a shows that both approximation and upper bound errors (”Approximation 2D” and ”Upper bound 2D”) decrease exponentially fast with the number of basis points. Moreover, they decrease with the same speed. This is an important property because in practice, we do not want to measure the true approximation error as it is very costly, but Fig. (a)a shows that the decrease of the upper bound error can be used as a convergence criterion for Algorithm 1.
We also tested our algorithm for 3d conductivity region, considering scalp conductivity as variable. The ”Upper bound 3D” in Fig. (a)a represents the decreasing of upper bound error with the number of basis points. We can notice that it is still exponentially fast, even if the absolute value of a slope is smaller.
Remark. In the case of symmetric BEM, head matrix is not homogeneous with respect to conductivities (even if the lead field is). So, considering scalp conductivity as unknown does increase the complexity of the projection problem (10).
We also evaluated the upper bound error on the points which did not belong to the sampling we used in algorithm, and was not significantly bigger then the one on sample points. This shows the good continuity properties of our method, mentioned previously.
IvC Approximation time
The time required for computing 625 exact lead field matrices is 4 hours, which results in 23 seconds per matrix (2 physical cores @ 2.60GHz, 16Gb RAM).
Precomputing of all required matrices on 20 basis points for Algorithm 2 is achieved in 14 minutes.
Once basis points are computed it takes 58 seconds to approximate 625 lead fields, which results in 0.09 second per matrix.
As mentioned in previous section, the approximation time does not depend on the complexity of the head model as it only requires to solve a leastsquares problem on the number of basis points. For more complex head models, the precomputing time will increase, but the approximation time will remain equal to 0.09 sec for 20 basis points.
IvD Conductivity estimation with simulated sources
We now examine convergence properties in a realistic application of conductivity estimation.
Using the same head model as above, we simulate a dataset which corresponds to a single dipole source for a reference conductivity (ground truth). Then, we use a simple dipole fitting approach to solve the inverse problem. For each of 625 conductivity samples in the domain, we approximate the lead field matrix and choose a source as a column of this matrix which fits measurement the best in terms of Euclidean norm. Let represent the th column of a matrix M. For each conductivity the data fitting error is defined as follows:
where denotes th source’s lead field, and  its optimal amplitude. represents the best fitting error of a single dipole corresponding to the measurement and conductivity .
Computing for each of 625 conductivity samples, we obtain a data fitting error map on the conductivity domain. We estimated the conductivity as one which minimizes this data fitting error.
For our simulation, we took the source whose lead field had the slowest convergence of the approximation error, i. e. the ”worst” source in terms of lead field approximation. Using the exact lead fields, we obtain a data fitting error map (Fig. 4) whose minimum lies at the simulated conductivity point (1.25 for brain and 0.05 for skull).
The goal here is to show that approximated lead fields are of sufficient quality for conductivity estimation. To do so, we compute lead field approximation for different number of basis points with our method. It is visible from Fig. (a)a that the ”shape” of the data fitting error map converges very fast to the exact one. 10 basis points are sufficient to correctly estimate the simulated conductivity.
IvE Conductivity estimation with real data
We use the same subject (head model) and conductivity domain of interest as in the previous section. Real EEG data was taken from a median nerve stimulation experiment. The right median nerve was percutaneously stimulated using monophasic squarewave impulses with a duration of 0.3 ms at 2.8 Hz. We filtered and averaged the data, and removed heartbeats and eye movements artifacts using Brainstorm software [21]. For more experiment and preprocessing details see [20].
We are interested in the N20P20 somatosensory averaged evoked potentials originating from Brodmann’s area 3b situated in the posterior bank of the Rolandic fissure [23, p. 1076]. We can see a remarkable activity peaking around 20 ms which has a dipolar topography on the left hemisphere (Fig. 5).
We analyzed 5 time samples () of a time window from 0.0185s to 0.0205s. The data fitting map is computed as follows:
The data fitting error, computed using the exact lead fields, is shown in Fig. 6. It is normalized to its minimum value, i.e. . Many factors contribute to the fitting error: additive noise, wrong conductivity model, simplified head and source models, inverse problem assumption, etc. The normalized data fitting error map can be interpreted as follows. Considering the error resulting from all these factors, except conductivity, to be equal to 1, the normalized error shows the relative impact of ”incorrect” conductivity estimation compared to other factors of error.
Because of the different sources of noise, the impact of the brain conductivity becomes less important and we can not significantly find its optimal value. But, for each fixed brain conductivity, we can define the skull conductivity which minimizes the error: this skull conductivity value lies between 0.01 S/m and 0.02 S/m. We can also notice that conductivity contributes up to 10% relatively to the other sources of error, which shows the importance of the correct estimation of this parameter.
As can be seen in Fig. (b)b, it is enough to use 11 basis points are sufficient to obtain an error map similar to the one obtained using exact lead fields.
V Conclusion
In this work, we introduced a method for fast approximation of the EEG forward problem solution for different conductivities. Computing the exact lead field for many conductivity configurations is time consuming, but our approach only needs to compute a relatively small number of exact solutions. All other lead fields can be approximated in a fast way with controlled approximation error. This accelerated computing time allows to explore conductivity space, which is crucial for noninvasive head tissue conductivity estimation from EEG data.
Our method provides both the lead field approximation with given basis points as well as the way to choose these points. It is done in a way to guarantee the monotonic convergence of approximation error. Moreover, we showed that the complexity of this approximation does not depend on the number of vertices in a head model.
Besides the theoretical properties of our algorithm, we studied its empirical performance. Realistic simulations showed the exponential decrease of approximation error with respect to the number of basis points. We also showed the usefulness of the method in realistic context of EEG inverse problem with both simulated and real EEG data. As expected, a relatively small number of precomputed basis points provide results which are similar but remarkably faster to using exact matrices.
The main motivation of this work was to propose a tool which would boost computing of lead fields, which is necessary for simultaneous estimation of brain sources’ activities and head tissues conductivity. We presented a simple approach of solving this inverse problem, based on a single dipole data fitting, to show that our method can be efficiently used for this kind of problems. We plan to continue studying the EEG inverse problem methods which consider conductivity as unknown but with more complex source models (e.g. multiple dipoles), with more unknown conductivities, e.g. composite skull structure [9], and different conductivity space sampling approaches, e.g. using gradient descent instead of uniformly sampled domain. We believe that our method will significantly improve the practical aspects of such studies.
Acknowledgment
This work was supported be ANR grant VIBRATIONS (ANR13PRTS0011) and has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (ERC Advanced Grant agreement No 694665 : CoBCoM  Computational Brain Connectivity Mapping).
References
 [1] S. Rush and D. A. Driscoll, “Current Distribution in the Brain From Surface Electrodes,” Anesthesia and Analgesia, vol. 47, pp. 717–723, 1968.
 [2] S. I. Gonçalves, J. C. De Munck, J. P. Verbunt, F. Bijma, R. M. Heethaar, and F. Lopes da Silva, “In vivo measurement of the brain and skull resistivities using an EITbased method and realistic models for the head,” IEEE Transactions on Biomedical Engineering, vol. 50, no. 6, pp. 754–767, 2003.
 [3] R. Hoekema, G. H. Wieneke, F. S. S. Leijten, C. W. M. Van Veelen, P. C. Van Rijen, G. J. M. Huiskamp, J. Ansems, and A. C. Van Huffelen, “Measurement of the conductivity of skull, temporarily removed during epilepsy surgery,” Brain Topography, 2003.
 [4] S. Vallaghé and M. Clerc, “A global sensitivity analysis of three and fourlayer EEG conductivity models,” IEEE Transactions on Biomedical Engineering, 2008.
 [5] Z. Akalin Acar and S. Makeig, “Effects of forward model errors on EEG source localization,” Brain Topography, 2013.
 [6] G. Wang and D. Ren, “Effect of braintoskull conductivity ratio on EEG source localization accuracy,” BioMed Research International, 2013.
 [7] R. Van Uitert, C. Johnson, and L. Zhukov, “Influence of head tissue conductivity in forward and inverse magnetoencephalographic simulations using realistic head models,” IEEE Transactions on Biomedical Engineering, 2004.
 [8] R. Pohlmeier, H. Buchner, G. Knoll, a. Rienäcker, R. Beckmann, and J. Pesch, “The influence of skullconductivity misspecification on inverse source localization in realistically shaped finite element head models.” Brain topography, 1997.
 [9] M. Dannhauer, B. Lanfer, C. H. Wolters, and T. R. Knösche, “Modeling of the human skull in EEG source analysis,” Human Brain Mapping, vol. 32, no. 9, pp. 1383–1399, sep 2011. [Online]. Available: http://doi.wiley.com/10.1002/hbm.21114
 [10] F. Costa, H. Batatia, T. Oberlin, and J.Y. Tourneret, “Skull Conductivity Estimation for EEG Source Localization,” IEEE Signal Processing Letters, vol. 24, no. 4, pp. 422–426, apr 2017. [Online]. Available: http://ieeexplore.ieee.org/document/7855647/
 [11] Z. Akalin Acar, C. E. Acar, and S. Makeig, “Simultaneous head tissue conductivity and EEG source location estimation,” NeuroImage, vol. 124, pp. 168–180, jan 2016.
 [12] S. Lew, C. H. Wolters, A. Anwander, S. Makeig, and R. S. MacLeod, “Improved EEG source analysis using lowresolution conductivity estimation in a fourcompartment finite element head model,” Human Brain Mapping, 2009.
 [13] S. Vallaghé, M. Clerc, and J. M. Badier, “In vivo conductivity estimation using somatosensory evoked potentials and cortical constraint on the source,” in 2007 4th IEEE International Symposium on Biomedical Imaging: From Nano to Macro  Proceedings, 2007.

[14]
J. S. Hesthaven, G. Rozza, and B. Stamm,
Certified Reduced Basis Methods for Parametrized Partial Differential Equations
, ser. SpringerBriefs in Mathematics. Cham: Springer International Publishing, 2016. [Online]. Available: http://link.springer.com/10.1007/9783319224701  [15] M. Hämäläinen, R. Hari, R. J. Ilmoniemi, J. Knuutila, and O. V. Lounasmaa, “Magnetoencephalography theory, instrumentation, and applications to noninvasive studies of the working human brain,” Reviews of Modern Physics, 1993.
 [16] J. Kybic, M. Clerc, T. Abboud, O. Faugeras, R. Keriven, and T. Papadopoulo, “A common formalism for the Integral formulations of the forward EEG problem,” IEEE Transactions on Medical Imaging, vol. 24, no. 1, pp. 12–28, jan 2005. [Online]. Available: http://ieeexplore.ieee.org/document/1375158/
 [17] C. H. Wolters, L. Grasedyck, and W. Hackbusch, “Efficient computation of lead field bases and influence matrix for the FEMbased EEG and MEG inverse problem,” Inverse Problems, 2004.
 [18] S. Vallaghé and T. Papadopoulo, “A Trilinear Immersed Finite Element Method for Solving the Electroencephalography Forward Problem,” SIAM Journal on Scientific Computing, vol. 32, no. 4, pp. 2379–2394, jan 2010. [Online]. Available: http://epubs.siam.org/doi/10.1137/09075038X
 [19] J. Sarvas, “Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem,” Physics in Medicine and Biology, 1987.
 [20] F. Tadel, Y. Haruta, E.i. Okumura, and T. Asakawa, “Yokogawa/KIT tutorial: Median nerve stimulation,” 2016. [Online]. Available: http://neuroimage.usc.edu/brainstorm/Tutorials/Yokogawa
 [21] F. Tadel, S. Baillet, J. C. Mosher, D. Pantazis, and R. M. Leahy, “Brainstorm: A userfriendly application for MEG/EEG analysis,” Computational Intelligence and Neuroscience, vol. 2011, pp. 1–13, 2011. [Online]. Available: http://www.hindawi.com/journals/cin/2011/879716/
 [22] A. Gramfort, T. Papadopoulo, E. Olivi, and M. Clerc, “OpenMEEG: opensource software for quasistatic bioelectromagnetics,” BioMedical Engineering OnLine, vol. 9, no. 1, p. 45, 2010. [Online]. Available: http://biomedicalengineeringonline.biomedcentral.com/articles/10.1186/1475925X945
 [23] E. Niedermeyer and F. Da Silva, Electroencephalography: Basic Principles, Clinical Applications, and Related Fields, Fifth Edition, 5th ed. Lippincott Williams and Wilkins, 2004.
Comments
There are no comments yet.