The estimation of the focal brain activity from electroencephalography (EEG) data is an ill-posed inverse problem, and the solution depends strongly on the accuracy of the head model, particularly the used skull conductivity value [1, 2]. The skull conductivity value can be calibrated by using well defined somatosensory evoked potentials / fields in combination with EEG , combined EEG/MEG  or by using electrical impedance tomography (EIT) -based techniques . However, these methods are sometimes unreliable, computationally effortful and may require auxiliary imaging tools.
. With BAE, the natural variability of the skull conductivity was modelled in a probability distribution, and subsequently its effects on the EEG data were taken statistically into account in the observation model (likelihood) with the help of an additive modelling error term. This allowed us to use an approximate skull conductivity in the EEG source imaging, i.e. to neglect the calibration step completely, and still get feasible source configuration estimates.
In this paper, we demonstrate with simulations that in addition to the source configuration the Bayesian uncertainty modelling also allows simultaneous estimation of the unknown skull conductivity. This estimation can be carried out through a low order decomposition of the modelling error term [8, 9].
We envision that the proposed approach could be used as such, or alternatively as a calibration step to determine the skull conductivity of a patient. The benefit compared to current injection -based techniques is that it relies solely on EEG data and thus avoids the requirement for a brain rest period before the actual EEG examination. In addition, this approach does not add computational effort (to the regular source imaging) since it does not involve non-linear formulations or unstable matrix inversions [3, 10]; instead, all the required pre-computations can be performed off-line before any experiments take place.
2.1 Bayesian framework with linear forward model
The computational domain is denoted with and its electric conductivity with where . The numerical observation model is
where are the measurements, is the number of measurements, is the (2-dimensional) lead field matrix that depends on , is the distributed dipole source configuration and is the measurement noise. Note that here the model assumes (unrealistically) that the accurate values of electric conductivities are known.
In the Bayesian framework, the inverse solution is the posterior density
where is the likelihood and the prior.
For the model (1), the likelihood can be written as
2.2 Bayesian approximation error (BAE) approach
In BAE, the accurate lead field matrix is replaced with an approximate (standard) lead field matrix in which fixed electric conductivities are used. Now, the observation model (1) can be written as
is the induced approximation error, .
In order to approximate the unknown skull conductivity, we need to obtain an estimate about the approximation error. Similarly as in [9, 8], we express as a linear combination of basis functions. Numerically, the estimation of
is tractable if most of its variance is captured in only few of the decomposition terms[9, 8].
A feasible set of basis functions that assigns most of the variance of
in the first few terms can be chosen based on the eigenvalue decomposing of the modelling error covariance. In particular, can be decomposed according to
where are the eigenvalues and
are the eigenvectors of. Based on this decomposition, we have that .
The approximation error can be expressed as
where and , with and . The coefficients and are given by the inner products and .
Thus, the observation model (4) can be rewritten as
where , contains the first eigenvectors and . Note that where .
For the simultaneous estimation of and , we have to construct the posterior model . To obtain a computationally efficient solution, we make the technical approximation that are mutually Gaussian and uncorrelated [9, 8]. Then, we obtain the approximate likelihood
where . Based on the eigenvalue decomposition and the Gaussian approximations, we have that the coefficients
are normally distributed with zero mean and covariance(where ). Thus, the posterior density becomes
where the Cholesky factors are and .
3.1 Head models
A 2-dimensional finite element (FE) -based three compartment (scalp, skull and brain) head model with 32 electrodes was used in the demonstrations (see, Figure 1). For the forward simulations and inverse estimations two different meshes with 2518 and 1780 nodes were constructed, respectively, with source spaces covering the gray matter of the brain.
In particular, we created head models with skull conductivity samples
drawn from a Gaussian distributionwith mean
(in S/m) and standard deviation 0.0013 (see, Figure1). We refer to these head models as sample head models. This skull conductivity distribution coincides with the literature values that range between 0.0041 and 0.033 [11, 12, 4]. The electric conductivity of the scalp was 0.43  and the brain 0.33 . We also created a standard head model with skull conductivity .
3.2 Computation of the approximation error statistics
The approximation error (5) depends both on the skull conductivity and the dipole source location. Feasible estimates for the skull conductivity can be expected if most of the variance of is related to the skull uncertainty. Therefore, we isolated the effect that the skull has on the observations. To do that, we estimated the error statistics for each dipole location separately. In the following text, the subscript denotes a specific dipole location.
The samples for the error statistics at location were created by first choosing randomly one of the sample head models and evaluating both the sample model and the standard model with the same single dipole source , i.e.
A set of single radial dipole samples with varying amplitudes were drawn from a Gaussian distribution with mean 1 and standard deviation 0.01, and 200 different leadfield matrices were used to estimate the error statistics at location . Superscript where and . The mean and error covariance at location were
Samples and were evaluated according to and , where and were matrices with columns the first and remaining eigenvectors of respectively.
3.3 Simultaneous skull conductivity and source estimation
In this work, we compute the maximum a posteriori (MAP) estimates of the posterior (10). Since we study only single dipole source cases, we can use a modified single dipole scan algorithm and solve
where is the total number of source locations, are the columns of the (standard) leadfield matrix corresponding to location , is the dipole reconstructed at location , are the basis functions obtained from the eigenvalue decomposition of the error covariance , are the reconstructed coefficients for these basis functions, and where parameters are the eigenvalues of the covariance (12). The number of parameters, as well as, the number of eigenvalues (used in the estimates) was defined based on the minimum value that satisfied . In most cases, either one or two parameters were needed.
The solution is the pair
where is the location that gives the lowest value for the functional (14).
Given the dipole location , the approximation error (5) can be written as
The skull conductivity can now be estimated by maximizing
where and are mutually independent.
From equation (15), we have that where the subscript is used to clarify that this is the probability density of , and the expression formally comes from the inversion of the lead field matrix difference and multiplication with . So, . Even though an have been modelled as Gaussian distributions, the non-linearity with respect to and often the infeasibility of estimating can prohibit a straightforward estimation. Instead we can use .
, the joint distribution ofand can be approximated as Gaussian, and therefore the MAP estimate (16) gives
where is the solution of problem (14), is the mean value of the postulated skull conductivity distribution, is the cross-covariance between and (estimated using samples and as described in section 3.2), and where are the first eigenvalues of the sampled-based error covariance (12).
4 Results, discussion and future work
To demonstrate the simultaneous dipole and skull conductivity reconstruction, we used two test cases with different skull conductivities. The EEG data was computed using the accurate
model that had either skull conductivity 0.0055 or 0.011 S/m, a single radial dipole source, and signal to noise ratio 30 dB. In Figure2, we show first the test set-up, then the dipole-scan reconstruction with the standard lead field that has skull conductivity 0.0085 S/m, and finally the proposed Bayesian error modelling result. To ease comparisons, we calculated the Euclidean distance (ED) (in milli meters) between the actual and reconstructed source. It can be seen that the proposed uncertainty modelling improves the source reconstructions when compared to the solution of the standard model. Moreover, the estimates for the skull conductivity are close to the true skull conductivities.
The feasibility of the skull conductivity estimates depends on how well the first few eigenvectors of describe the conductivity related modelling errors. In the future, we aim to isolate the conductivity related errors from other error sources (e.g. discretization) in order to improve the skull conductivity estimates.
In our analysis the head geometry was considered fully known; however, it can also be taken into account as an uncertainty . We also note that in this work we estimated a single conductivity value for the whole skull. As a next step, we plan to use the proposed approach in conjugation with other optimization techniques [15, 16] to produce more detailed skull conductivity maps.
We have demonstrated that with the help of the Bayesian uncertainty modeling the EEG source imaging problem could be reformulated in such a way that it allows to estimate, simultaneously with the source configuration, a low-order approximation for the modelling error. With the help of this modelling error estimate, the underlying skull conductivity could be reconstructed.
This project has received funding from the ATTRACT project funded by the EC under Grant Agreement 777222 and from the Academy of Finland post-doctoral program (project no. 316542).
Conflict of interest
The authors declare that they have no conflicts of interest.
-  B. Vanrumste, G.V. Hoey, R.V. de Walle, M. D’Have, I. Lemahieu, P. Boon, Med. Biol. Eng. Comput. 38(5), 528 (2000).
-  M. Dannhauer, B. Lanfer, C. Wolters, T. Knösche, Hum. Brain Mapp. 32, 1383 (2011)
-  S. Lew, C.H. Wolters, A. Anwander, S. Makeig, R.S. MacLeod, Hum. Brain Mapp. 30(9), 2862 (2009).
-  U. Aydin, J. Vorwerk, P. Küpper, M. Heers, H. Kugel, A. Galka, L. Hamid, J. Wellmer, C. Kellinghaus, S. Rampp, C.H. Wolters, PLoS ONE 9, e93154 (2014)
-  S.I. Goncalves, J.C. de Munck, J.P. Verbunt, R.M. Heethaar, F.H. da Silva, IEEE Trans. Biomed. Eng. 50(9), 1124 (2003)
-  V. Rimpiläinen, A. Koulouri, F. Lucka, J.P. Kaipio, C.H. Wolters, IFMBE Proc. 65, 892 (2018)
-  V. Rimpiläinen, A. Koulouri, F. Lucka, J.P. Kaipio, C.H. Wolters, NeuroImage 188, 252 (2019).
-  J. Kaipio, V. Kolehmainen, in Bayesian Theory and Applications, ed. by P. Damien, N. Polson, D. Stephens (Oxford University Press, 2013)
-  A. Nissinen, V. Kolehmainen, J.P. Kaipio, Int. J. Uncertain. Quantif. 1(3), 203 (2011)
-  C. Papageorgakis, Patient specific conductivity models: characterization of the skull bones. Ph.D. thesis (2017)
-  R. Hoekema, G.H. Wieneke, C.W. van Veelen, P.C. van Rijen, G.J. Huiskamp, J. Ansems, A.C. van Huffelen, Brain Topogr. 16(1), 29 (2003)
-  S. Homma, T. Musha, Y. Nakajima, Y. Okomoto, S. Blom, R. Flink, K.E. Hagbarth, Neurosci. Res. 22(1), 51 (1995)
-  C. Ramon, P.H. Schimpf, J. Haueisen, Biomed. Eng. Online 5(10) (2006)
-  A. Koulouri, V. Rimpiläinen, M. Brookes, J.P. Kaipio, Appl. Num. Math. 106, 24 (2016)
-  Z.A. Acar, C.E. Acar, S. Makeig, NeuroImage 124, 168 (2016).
-  F. Costa, H. Batatia, T. Oberlin, J.Y. Tourneret, arXiv: 1609.06874 (2016)