Forward and Inverse Effects of the Complete Electrode Model in Neonatal EEG

!! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Abstract ! This paper investigates finite element method (FEM) based modeling in the context of neonatal electroencephalography (EEG). In particular, the focus lies on electrode boundary conditions. We compare the complete electrode model (CEM) to the point electrode model (PEM), which is the current standard in EEG. In the CEM, the voltage experienced by an electrode is modeled more realistically as the integral average of the potential distribution over its contact surface, whereas the PEM relies on a point value. Consequently, the CEM takes into account the sub-electrode shunting currents which are absent in the PEM. In this study, we aim to find out how the electrode voltage predicted by these two models differ, if standard size electrodes are attached to a head of a neonate. Additionally, we study voltages and voltage variation on electrode surfaces with two source locations: (A) next to the T8 electrode and (B) directly under the Fz electrode and the frontal fontanel. A realistic model of a neonatal head including a skull with fontanels and sutures is used. Based on the results, the forward simulation differences between CEM and PEM are in general small, but significant outliers can occur in the vicinity of the electrodes. The CEM 1(Corresponding author, Email: Sampsa.Pursiainen@tut.fi). Department of Mathematics, Tampere University of Technology, Finland. 2Newborn Medicine in the Boston Children’s Hospital, MA, USA. 3Department of Engineering, Olivet Nazarene University, IL, USA. 4Institute for Biomagnetism and Biosignalanalysis, University of Münster, Germany.


Introduction
This paper explores the mathematical complete electrode model (CEM) (Cheng et al., 1989;Somersalo et al., 1992;Ollikainen et al., 2000;Vallaghé et al., 2009;Pursiainen, 2012;Vauhkonen, 1997) in neonatal electroencephalography (EEG) forward simulation (Lew et al., 2013;Niedermeyer and da Silva, 2004;Gargiulo et al., 2015;Baillet et al., 2001;de Munck et al., 2012).Our objective is to evaluate the differences between the CEM and the standard point electrode model (PEM) for a head of a three days old neonate for which the electrode diameter is unusually large compared to that of the head.The CEM consists of a set of boundary conditions originally developed for impedance tomography (Cheney et al., 1999) in which the electrodes are utilized to inject currents and measure voltage data simultaneously, thus necessitating accurate control of impedances and surface currents to minimize forward simulation errors.
In EEG, the net currents flowing through the electrodes are zero.Sub-electrode surface currents due to potential variation, however, exist.These shunting effects consume a small part of the electric field energy, leading to lower voltage amplitudes than what is predicted by a mathematical model not including the skinelectrode interaction.The shunting is the more pronounced the larger the relative electrode size (angular coverage) or the lower the contact impedance is (Ollikainen et al., 2000).The first aspect of these makes the neonatal EEG interesting as a target of investigation.The second one is important from the viewpoint of choosing an optimal measurement impedance: the general rule of thumb is that decreasing the impedance improves the signal amplitude, because of better power transfer between the skin and the electrode, until around 100-500 Ohm is reached.Below 100 Ohm, shunting increases significantly altering the electrode voltage and the surface potential distribution (Duffy et al., 2012;American Clinical Neurophysiology Society, 2006).
The CEM accounts for more accurate modeling of electrode shunting, and thereby, in principle, extends the interval of applicable impedances beyond that of the classical PEM.In the CEM, the voltage experienced by an electrode is the integral average of the potential distribution below its contact surface, whereas the PEM relies on a single mesh node.This study aims to find out how the electrode voltage predicted by these two models differ, when 10 mm diameter electrodes are used in combination with a newborn head.
Additionally, the voltage values and variation on electrode surfaces are investigated utilizing two source locations: (A) next to the T8 (5-th) electrode and (B) directly under the Fz electrode and the frontal fontanel.Because of the significant differences between a newborn and a normal adult skull, this study utilizes a realistic neonatal head model including a skull with fontanels and sutures (Lew et al., 2013, Darbas et al. 2016).As the primary forward simulation method, we use the finite element method (FEM) (de Munck et al., 2012;Lew et al., 2009;Braess, 2001), which is a flexible alternative to the boundary element method (Mosher et al., 1999;Acar and Makeig, 2010;Gramfort et al., 2011;Stenroos and Sarvas, 2012) as well as to finite difference and volume methods (Montes-Restrepo et al., 2014;Wendel et al., 2008;Vatta et al., 2009;Cook and Koles, 2006).Particularly important for our study is that a FEM mesh can be simultaneously adapted to both boundary and interior structures of the computational domain.Our CEM implementation is mathematically rigorous: The lead field matrix describing the forward simulation is constructed directly based on the weak form of the electric potential (Poisson) equation.Moreover, the source current distribution is modeled via the Whitney (Raviart-Thomas) type vector basis functions with divergence satisfying square integrability condition (Bauer et al., 2015;Pursiainen, 2012;Pursiainen et al., 2011;Tanzer et al., 2005).
This study includes numerical experiments in which a wide range of electrode impedances from 0.1 Ohm to 12 kOhm are investigated.The differences of the electrode voltages predicted by CEM and PEM are explored utilizing a descriptive statistical approach in which the relative difference and magnitude measures (RDM and MAG) are analyzed via box-plots for five different source eccentricities (relative norms).Additionally, voltages and voltage variation on electrode surfaces with two different electrode impedance values are visualized.The inverse aspect of the forward simulation differences is explored through the minimum current estimation (MCE) technique (Calvetti et al., 2009;Lucka et al., 2012;Uutela et al., 1999).
Based on our results, the CEM can be considered as an integral part of the head model.The median level of the RDM and MAG between CEM and PEM was found to be small, around 0.7 and 0.5 %, regarding the impedances above the commonly recommended lower limit 100 Ohm.However, the dependence of these differences on the source position wasfound to be significant: the spread of the results extended considerably towards the electrodes with outliers up to 23 % and 36 % occurring in their vicinity.

Forward Model
Let e  ,  =1, 2, …, L ( L = the number of electrodes) denote the skin-electrode contact surfaces on the boundary ∂Ω of the head Ω .Given the primary current density  j p and a conductivity tensor σ within Ω , the resulting electric potential field u within Ω can be obtained as a solution of the quasi-static electric potential equation with the following CEM boundary conditions (Cheng et al., 1989): σ ∇ (4) The first one of these requires that the normal current σ ∇u ⋅  n on ∂Ω can flow out of or into the domain only through electrodes.The second one ensures that the net current flowing through each electrode is zero, and the third one defines the potential jump on the skin-electrode contact interface.The potential jump is proportional to the normal current density and to the pointwise effective skin-electrode contact impedance  Z  with the unit Ohm m 2 .The voltage of the k-th electrode is denoted by U  .For simplicity, we assume that  Z  is of the form  Z  = Z  A  , where Z  is the average impedance (Ohm) of the contact surface area A  (m 2 ).Under this assumption, the ℓ-th electrode voltage is given by the integral mean ∫ dS.

Weak Form
The weak form (Somersalo et al., 1992;Vauhkonen, 1997) of ( 1)-( 4), i.e., the CEM, can be written as follows (Appendix): (5) Here, dV and dS denote differentials with respect to volume and surface, respectively.If the divergence of  j p is square integrable, i.e., if it holds that ( ) holds for u,v,uv, and U l , and consequently, the CEM weak form tends to PEM one, i.e., the finite element (FE) mesh T, the potential and primary current density can be approximated as the following finite sums u T = within the subspace of the basis functions is given by Matrix A is of the form where additionally a i', j = δ i', j with δ i',i' =1 and δ i', j = 0 for i' ≠ j (Kronecker delta) for a single nodal basis function ψ i' attaining its maximum on the part of the boundary not covered by the electrodes.The Electrode voltages y = U 1 , U 2 , …, U L ( ) predicted by ( 7) can be obtained via y = R v, in which multiplication by the real L-by-L matrix R, r , =1-1 / L for  =1, 2,…, L and r i, = −1 / L for i ≠ , fixes the sum of the electrode voltages to zero.Consequently, the potential data predicted by the forward model can be simulated as given by y = Lx, where the lead-field matrix is of the following form In order to determine the lead-field matrix L efficiently, the matrix following from (6).

Synthetic Dipoles
The primary current and electric potential fields were modeled utilizing a Whitney (Raviart-Thomas) and linear Lagrangian (nodal) function basis, respectively.This combination yields a simple model of synthetic dipoles (Bauer et al., 2015;Pursiainen et al., 2011;Pursiainen, 2012;Tanzer et al., 2005) with dipole moment and position determined by where  r P j and  r P i are the position vectors of mesh nodes P j and P i on the opposite sides of the common face F in an adjacent pair of tetrahedra (Figure 1).The synthetic dipole moment is given by the integral dV and the position is the midpoint of the nodes P j and P i for which with s ψ ,P = 1, if ψ corresponds to P, 0, otherwise.

Minimum Current Estimation
The inverse aspect of modeling differences was explored via the MCE approach (Uutela et al., 1999) that where γ > 0 is a regularization parameter.The minimizer was estimated using the following iterative alternating sequential (IAS) algorithm (Calvetti et al., 2009): 2. Find x i+1 as the least-squares solution of the linear system where ) , and set i → i +1.
3. If i is less than the desired number of iterations, repeat the step number two.

4.
Otherwise, set x = x i as the final estimate.
The result is a minimally supported current distribution x (Appendix) which is here interpreted as a single dipole.Based on x, estimates of dipole location and direction were obtained as the weighted average

Numerical Experiments
Our numerical experiments utilized three variations (I)-(III) of a neonate head model which was created based on the segmentation of T1 weighted MPRAGE images of a 3-day-old healthy infant (Lew et al. ,2013).The segmentation was conducted manually in FreeView, a volume and surface visualization tool within the FreeSurfer software application (Dale et al., 1999), resulting in segmented skull, scalp, cerebrospinal fluid, grey and white matter, and anterior fontanel compartments.Other fontanels including the posterior, sphenoidal, and mastoid fontanels were segmented based on the knowledge of skull geometry was generated for the neonate segmentation (Figure 2).
Model (I) was the gold standard with the following five tissue compartments and values of conductivity (Table 1): the brain (grey and white matter) 0.33 S/m (Haueisen, 1996), cerebrospinal fluid (CSF) 1.79 S/m (Baumann et al., 1997), bone 0.04 S/m and fontanel/suture stuctures (openings) within the skull 0.3 S/m (Lew et al., 2013), and skin 0.33 S/m (Haueisen, 1996).Model (II) was otherwise identical to (I) except that the conductivity of fontanel/suture (0.3 S/m), characteristic to a neonatal head, was replaced with that of bone (0.04 S/m), resulting in a closed skull.Models (I) and (III) differed with respect to the bone conductivity, which in (III) was given the value 0.1 S/m.The voltage data was obtained utilizing a cap of 74 circular electrodes with 10 mm diameter (10-20 system).Each electrode was formed as a set of surface triangles with center point within 5 mm distance from the electrode center.In the visualization of the results, the nearest point interpolation method between the original and a refined surface mesh was utilized.Fourteen different impedance values, equally distributed on a logarithmic scale, were used to cover a range of values between 0.10 Ohm and 12 kOhm average contact impedance (ACI).
The relative difference and magnitude (RDM and MAG) measures between two voltage predictions u 1 vs. u 2 were evaluated.RDM approximates topography changes and MAG the difference in magnitude between two forward solutions.The electrode/head model comparisons investigated are listed in Table 2. Comparison (1) is the primary one and ( 2  r w and  r as in ( 11) and ( 16), respectively, were analyzed via source localization analysis in which each one of the sources was reconstructed with the IAS algorithm.For each source, the exact voltage data were calculated utilizing the CEM and the reconstruction was computed via PEM forward simulation.The synthetic dipole used in generating the data was not included in the reconstruction procedure in order to avoid the inverse crime condition, i.e., overly optimistic data fit.
The number of IAS iteration steps was chosen to be 50 and the regularization parameter γ was given the value 1.4•10 −6 .The present choice of γ can be motivated through a correspondence to a hierarchical Bayesian probability model with a conditionally Gaussian prior and a gamma hyperprior (Calvetti et al., 2009;Lucka et al., 2012): Assuming that the measurements contain Gaussian white noise with very low standard deviation ν = 0.001 and setting the scaling (initial prior variance) and shape parameter of the gamma density to θ 0 =1 and β =1.5, then the resulting posterior probability density can be maximized via the IAS algorithm by setting γ = v 2 2 /θ 0 equal to 1.4• 10 −6 (Calvetti et al., 2009).
Additionally, electrode voltages and voltage variation u −U  underneath the electrodes  =1, 2, …, L were analyzed for two 10 nAm sources (A) and (B) placed next to the T8 (5-th) electrode and directly under the Fz electrode and the frontal fontanel, respectively.The voltages corresponding to (A) and (B) were produced using the following electrode/head model and impedance combinations: PEM (I), CEM (I) 2.0 kOhm, CEM (I) 0.1 Ohm, and CEM (II) 2.0 kOhm.

Results
The results have been included in Figures 3-8.Of these, Figure 3  In comparison (1), the median RDM and MAG vary from 1.5 to 0.7 % (RDM) and -3.5 % to -0.5 % (MAG) when moving from 0.1 Ohm to 12 kOhm.The maximal (absolute) value decreases from 27 % to 23 % (RDM) and from 48 % to 36 % (MAG), respectively.Based on both RDM and MAG, the difference distributions stay virtually unchanged between 100 Ohm and 12 kOhm.The point of the steepest slope is located between 10 and 100 Ohm.In comparisons ( 1) and ( 3), the median and spread of both RDM and MAG coincide up to the accuracy of 0.2, whereas in (2) those appear to be on a lower level.
Figure 4 analyzes the results of comparison (1) with respect to the source depth for eccentricities.Each bar visualizes a distribution of 100 random sources.The absolute values of both RDM and MAG can be observed to grow along with the eccentricity.Close to the surface (eccentricity 0.98), the medians of RDM and MAG were found to vary between 2 -4 % and -1 --11 %, respectively.
PD, AD and ND calculated for comparison (1) were observed to increase in absolute value towards the surface of the head (Figure 5).At 0.98 eccentricity, the median PD was 2.1 -3.0 mm and the median AD 4 -6 degrees.As indicated by the negative median ND of -2.8 --1.8 mm, the sources were localized deeper than they actually were.ND and PD are close in magnitude, suggesting that the localization difference is mainly in the normal direction.Figure 6 suggests that large values in the absolute value of RDM/MAG correlate with those of PD/AD/ND in a non-systematic way.

Discussion
Our numerical experiments focused on forward simulation differences between the complete electrode model (CEM) and the point electrode model (PEM) boundary conditions with a realistic neonatal head model as the target domain.These experiments were motivated by the small size of the head compared to which the typical electrode diameter of 10 mm is large.Consequently, the differences can be expected to be larger than in a normal adult head.In order to investigate these changes, we investigated the relative difference and magnitude measures (RDM and MAG) for electrode impedances between extremely low 0.1 Ohm and comparably high 12 kOhm as well as for five different source eccentricities.Minimum current estimation (MCE) based source localization was utilized to enlighten the effect of forward simulation differences and also visual comparisons of the voltage patterns and sub-electrode potential distributions were presented.
The results suggest that the median level of RDM and MAG appears to be rather small with 0.5 -3.5 % for all impedances.The outliers were significant in the vicinity of the electrodes throughout the interval of tested impedances, extending in absolute value to 23 and 36 % for RDM and MAG, respectively, for the normally used range of over 100 Ohm.Hence, it is obvious that the CEM can improve the forward simulation accuracy in cortical regions.Above 100 kOhm the impedance value has a negligible effect on RDM and MAG.The exact value of impedance might thus be of minor importance in applications of experimental EEG with respect to the forward simulation accuracy.Based on the MAG, it is obvious that the PEM led to systematically lower amplitudes than the CEM due to the absence of the shunting effects.
As a whole, the differences between CEM and PEM (comparison (1)) were found to be similar, in terms of descriptive statistics, as those between the open (I) vs. closed (II) skull model (comparison (2) and ( 3)).In both cases large values of comparable magnitude were observed in the vicinity of the skull.In the former case, these outliers occur close to the electrodes, and in the latter one, close to the openings, i.e., the fontanels and sutures.Due to the difficulty of creating a point cloud without a bias in either of these regions, the eccentricity dependence was analyzed only for comparison (1), which was of the primary interest.The CEM-PEM differences were also observed to increase along with the conductivity of the skull tissue based on the results obtained with 0.04 and 0.1 S/m of (I) and (III).Hence, the CEM can be considered as an integral part of the outer head model akin to accurate skull modeling (Lew et al., 2013).
The MCE results suggest that the differences between the CEM and PEM can be reflected in terms of source localization results.An increase in RDM and MAG between CEM and PEM leads to slightly larger localization differences PD, AD and ND.With respect to impedance, marginal deterioration of the localization accuracy was observed below 100 Ohm.Above that value, the median PD was maximally 2.5 mm, which is close to the real extent of the neural sources (de Munck et al., 2012).The median ND was below zero close to the electrodes, suggesting that PEM based inverse estimates tend to localize sources 1-2 mm deeper than they actually are.This is obviously connected to the absence of the shunting effects in PEM, which is why the actual voltage magnitude is, in general, lower than what is predicted by the PEM.
The general level of AD was significant, which is at least partially due to the MCE inversion strategy being non-optimal with respect to finding the direction.However, the distribution of AD was similar to that of PD regarding the eccentricity and impedance, showing that also with respect to AD the accuracy of the MCE inversion decreases as the source modeling differences go up.MCE was chosen as the inverse approach, since it does not need the a priori knowledge of the number of underlying sources and can thus often better be used to recover brain activity in a real application context.
In order to assure the reliability of the results a realistic head model and an advanced source modeling technique were utilized in the numerical experiments.The accuracy of the applied Whitney type model has been recently shown to be slightly superior to the classical St. Venant (Schönen et al., 1994;Buchner et al., 1997;Toupin, 1965) and Partial Integration (Yan et al., 1991;Weinstein et al., 2000;Awada et al., 1997) approaches at the orientations and positions of the synthetic dipoles and between these two, slightly inferior to St. Venant, if arbitrary sources are in question (Bauer et al., 2015).Nodes connected to the surface of the brain were excluded from the forward simulation in order to eliminate outliers due to material parameter discontinuities (de Munck et al., 2012).The head model involves some uncertainty, since a neonatal head changes dramatically within months after the birth and especially changes in the skull compartment have significant impact on the EEG forward problem (Dannhauer et al., 2011, Lew et al., 2009).In particular, rapid ossification quickly decreases the conductivity of the skull, a higher value 0.1 S/m was tested in addition to 0.04 S/m (Lew et al., 2013).An increase in the conductivity of the skull was observed to slightly strengthen the differences between CEM and PEM.It is noteworthy that a recent study (Gargiulo et al., 2015) suggests than an even higher skull conductivity could be considered, which would lead to even higher errors.
Significant for clinical neonatal studies is the finding that the electrode and skull modeling errors might be comparable with each other.It has been shown that the neonatal unossified skull structure allows improved source localization when compared to an adult EEG because of the higher focality of dipolar EEG patterns due to higher skull conductivity (Odabaee et al., 2014).In order to achieve optimal source localization results, accurate electrode modeling might be necessary.An important aspect regarding clinical applications is also that the CEM can be implemented in a straightforward way via the approach presented in this paper.
It can be easily included to any EEG forward solver by associating an electrode with a set of elements instead of a single point and by forming the boundary condition matrices B and C that are not present in the PEM.Since it seems that the CEM can improve the modeling accuracy as compared to the PEM, the CEM is a well-motivated part of the forward solver.
To conclude our study, we propose that the CEM can improve EEG forward simulation accuracy.Yet, it is not a completely necessary feature in source modeling of standard EEG measurements, where the impedance is over 100 Ohm.In our future research, the CEM will be utilized to combined transcranial current stimulation (tCS) (Herrmann et al., 2013;Dannhauer et al., 2012;Eichelbaum et al., 2014;Breitling et al., 2016;Guler et al., 2016b,a) and EEG, in which the CEM can be used for both stimulation and measurement electrodes.As motivated by the current results, numerical tests could be performed over a range of skull sizes and parameter values.An important topic would be also to study the CEM utilizing experimental data, in order to validate the present numerical results in practice.A well-controlled phantom (Liehr and Haueisen, 2007;Wetterling et al., 2009), or possibly even more informative, an animal model study (Lau et al., 2016) should be one of the next steps in the further evaluation of the CEM.In (Lau et al., 2016), the important first steps have been carried out in characterizing errors in source reconstruction due to ignoring skull defects and in assessing the ability of an exact FE head model to eliminate such errors in an evaluation study with an anesthetized rabbit.6. Appendix

Integration by parts
Equation ( 1) multiplied with a test function v ∈ S and integrated by parts leads to The last term on the right-hand side can be written as follows: This can be expanded as The final form (24) substituted into (22) yields (5).
Model   q w are determined by the midpoint and directional unit vector of the line segment between P i and P j , respectively.(Bauer et al., 2015).
) and (3) are used as a reference in order to show the voltage differences obtained with the open (I) and closed (II) skull model as well as the effect of alternative skull conductivity (III).The results were visualized via box-plots showing the median, spread (IQR), i.e., the centermost (thickened) interval containing 50 % of the sample points, as well as the (narrow) interval between the minimal and maximal result.Further comparison between CEM (I) and PEM (I) concerned random samples of 100 random sources collected at five different source depths in the brain compartment with eccentricities (relative norms) 0.2, 0.4, 0.6, 0.8, and 0.98.In addition to evaluation of RDM and MAG, the position, angular (degree) and norm (depth) difference (PD, AD and ND) with PD r w ,  Figures 4 and 5, respectively, and the dependence between RDM/MAG and PD/AD/ND is illustrated in Figure 6.Finally, Figures 7 and 8 visualize the electrode voltages and sub-electrode potential variation for

Figure 7
Figure7shows voltage patterns for sources (A) and (B) as obtained with the modeling techniques PEM (I), CEM (I) 2.0 kOhm, CEM (I) 0.1 Ohm, and CEM (II) 2.0 kOhm.In the case of (A), PEM (I) produces an outlier at the T8 electrode whereas the other techniques result in mutually similar patterns.In the case of (B), CEM (II) 2.0 Ohm yields around 1-2 µV lower voltage amplitudes than the other techniques for the AFz and F2 (20-th and 69-th) electrode which are close to the source.Figure8illustrates the voltage variation beneath the electrodes for (A) and (B), showing that the incorrect approximation obtained with PEM (I) for (A) at the 5-th electrode appears to be due to the difference between the electrode center point SP was supported by the Centre of Excellence in Inverse Problems and the key project no.305055 of the Academy of Finland.CHW was supported by EU project ChildBrain (Marie Curie Innovative Training Networks, no.641652) and by the Priority Program 1665 of the Deutsche Forschungsgemeinschaft (DFG) (WO1425/5-1,5-2,7-1).

Figure 1 :
Figure 1: A schematic illustration of the synthetic dipole corresponding to the face F .The position  r w and

Figure 2 :
Figure2: The neonatal head model.On left: Axial view of the skull (blue) together with fontanel/suture (orange).On right: Sagittal view of the grey matter (green) and fontanel/suture (orange).

Figure 3 :
Figure 3: The distributions of RDM (top) and MAG (bottom) for the electrode/head model comparisons (1) CEM (I) vs. PEM (I), (2) CEM (I) vs. CEM (II) and (3) CEM (III) vs. PEM (III) covering the tested range of impedances from 0.1 Ohm to 12 kOhm.The horizontal axis has a logarithmic scaling and unit Ohm.The unit of the vertical axis is %.

Figure 6 :
Figure 6: PD/AD/ND plotted against RDM/MAG suggests that large values in the absolute value of RDM/MAG correlate with those of PD/AD/ND.Outliers in RDM/MAG tend to result in large deviations of PD/AD/ND.Each point cloud includes all the results obtained for a given impedance value and eccentricities (relative norms) 0.2, 0.4, 0.6, 0.8, and 0.98.

Figure 7 :
Figure7: Voltage variation for two 10 nAm sources (A) and (B).The most significant variation occurs on the electrodes nearby the source.The variation can be observed to diminish along with the impedance due to increasing shunting effects: for 0.1 Ohm impedance, it is close to zero.

Figure 8 :
Figure 8: The voltage data for sources (A) and (B) located next to the T8 (5-th) electrode and below the Fz electrode and the frontal fontanel, respectively.In the case of (A), PEM (I) produces an outlier at the T8 electrode.Also in the case of (B), PEM (I) yields significantly different values compared to other techniques for electrodes close to the source, e.g.AFz and F2 (20-th and 69-th) electrode.