|Home | About | Journals | Submit | Contact Us | Français|
Bioelectric source analysis in the human brain from scalp electroencephalography (EEG) signals is sensitive to geometry and conductivity properties of the different head tissues. We propose a low resolution conductivity estimation (LRCE) method using simulated annealing optimization on high resolution finite element models that individually optimizes a realistically-shaped four-layer volume conductor with regard to the brain and skull compartment conductivities. As input data, the method needs T1- and PD-weighted magnetic resonance images for an improved modeling of the skull and the cerebrospinal fluid compartment and evoked potential data with high signal-to-noise ratio (SNR). Our simulation studies showed that for EEG data with realistic SNR, the LRCE method was able to simultaneously reconstruct both the brain and the skull conductivity together with the underlying dipole source and provided an improved source analysis result. We have also demonstrated the feasibility and applicability of the new method to simultaneously estimate brain and skull conductivity and a somatosensory source from measured tactile somatosensory evoked potentials of a human subject. Our results show the viability of an approach that computes its own conductivity values and thus reduces the dependence on assigning values from the literature and likely produces a more robust estimate of current sources. Using the LRCE method, the individually optimized four-compartment volume conductor model can in a second step be used for the analysis of clinical or cognitive data acquired from the same subject.
The electroencephalographic inverse problem aims at reconstructing the underlying current distribution in the human brain using potential differences measured non-invasively from the head surface. A critical component of source reconstruction is the head volume conductor model used to reach an accurate solution of the associated forward problem, i.e., the simulation of the electroencephalogram (EEG) for a known current source in the brain. The volume conductor model contains both the geometry and the electrical conduction properties of the head tissues and the accuracy of both parameters has direct impact on the accuracy of the source analysis [Buchner et al., 1997; Huiskamp et al., 1999; Gençer and Acar, 2004; Ramon et al., 2004; Zhang et al., 2006,2006b,2008; Rullmann et al., 2008]. The practical challenges of creating patient specific models currently prohibit this degree of customization for each routine case of clinical source analysis, thus it is essential to identify the parameters that have the largest impact on solution accuracy and to attempt to customize them to the particular case.
Magnetic Resonance (MR) or Computed Tomography (CT) imaging provides the geometry information for the brain, the cerebrospinal fluid (CSF), the skull, and the scalp [Huiskamp et al., 1999; Ramon et al., 2004; Wolters et al., 2006; Zhang et al., 2006b, 2008]. MRI has the advantage of being a completely safe and noninvasive method for imaging the head, while CT provides better definition of hard tissues such as bone. However, CT is not justified for routine physiological studies in healthy human subjects. In this study we used a combination of T1-weighted MRI, which is well suited for the identification of soft tissues (scalp, brain) and proton-density (PD) weighted MRI, enabling the segmentation of the inner skull/outer CSF surface. This approach leads to an improved modeling of the CSF compartment and of the skull thickness over standard (T1) weighted MRI, important for a successful application of the proposed low resolution conductivity estimation (LRCE) method. The volume conductor model used in this study consisted of four individually and accurately shaped compartments, the scalp, skull, CSF, and brain.
Determining the second component of the head model, the conductivities of the tissues, does not have the support of a technology as capable as MRI or CT. First attempts to measure the conductivities of biological tissues were in vitro, often using samples taken from animals [Geddes and Baker, 1967]. The conductivity of human CSF at body temperature has been measured by [Baumann et al., 1997] to be 1.79 S/m (average over 7 subjects, ranging in age from 4.5 months to 70 years, with a standard deviation of less than 1.4% between subjects and for frequencies between 10 and 10,000Hz). Based on the very low standard deviation determined in this study, the conductivity of the CSF compartment will be fixed to 1.79S/m throughout this study. EEG measurements were furthermore shown to be sensitive to a correct modeling of this highly conducting compartment, which is located between the sources in the brain and the measurement electrodes on the scalp ([Huang et al., 1990; Ramon et al., 2004; Wolters et al., 2006; Wendel et al., 2008; Rullmann et al., 2008], see also discussion and further references in [Baumann et al., 1997]). In contrast, the electric conductivities of skull and brain tissues were shown to vary much stronger across individuals and within the same individual due to variations in age, disease state and environmental factors. [Latikka et al., 2001] investigated the conductivity of living intracranial tissues from nine patients under surgery. As the skull has considerably higher resistivity than the other head tissues—and thus could be expected to play an especially big role in the electric currents in the head—much attention has been focused on determining its conductivity. Rush and Driscoll measured impedances for a half-skull immersed in fluid [Rush and Driscoll, 1968, 1969] and since then the brain:skull conductivity ratio (in three-compartment head models) of 80 has been commonly used in bioelectric source analysis [Homma et al., 1995]. A similar ratio of 72 averaged over six subjects was reported recently using two different in vivo approaches [Gonçalves et al., 2003a], one method using the principles of electrical impedance tomography (EIT) and the other method based on an estimation through a combined analysis of the evoked somatosensory potentials/fields (SEP/SEF. However, those results remain controversial because other studies have reported the following ratios: 15 (based on in vitro and in vivo measurements) [Oostendorp et al., 2000], 18.7 ± 2.1 (based on in vivo experiments using intracranial electrical stimulation in two epilepsy patients) [Zhang et al., 2006], 23 (averaged value over nine subjects estimated from combined SEP/SEF data) [Baysal and Haueisen, 2004], 25 ± 7 (estimated from intra- and extra-cranial potential measurements) [Lai et al., 2005], and 42 (averaged over six subjects using EIT measurements) [Gonçalves et al., 2003b]. At this point, there is little hope of a resolution of these large discrepancies, some of which may originate in inter-patient differences or natural variations over time (see, e.g. [Haueisen, 1996; Goncalves et al., 2003b]), some might result from ignoring the high conductivity of the CSF since most of the above studies used three-compartment (scalp, skull, brain) head models or from ignoring the influence of realistic geometrical shape when using spherical head models, so that we propose a four-compartment realistically-shaped head modeling approach that seeks to resolve variation for each individual case by making skull and brain conductivity an additional parameter to be solved.
The growing body of evidence suggesting that the quality and fidelity of the volume conductor model of the head plays a key role in solution accuracy [Cuffin, 1996; Huiskamp et al., 1999; Ramon et al., 2004; Rullmann et al., 2008] also drives the choice of numerical methods. There is a wide range of approaches including multi-layer sphere models [de Munck and Peters, 1993], the boundary element method (BEM) [Sarvas, 1987; Hämäläinen and Sarvas, 1989; de Munck, 1992; Fuchs et al., 1998; Huiskamp et al., 1999; Kybic et al., 2005], the finite difference method (FDM) [Hallez et al., 2005] and the finite element method (FEM) [Bertrand et al., 1991; Haueisen, 1996; Marin et al., 1998; Weinstein et al., 2000; Ramon et al., 2004; Wolters et al., 2006; Zhang et al., 2006, 2006b, 2008]. The FEM offers the most flexibility in assigning both accurate geometry and detailed conductivity attributes to the model at the cost of both creating and computing on the resulting geometric model. The use of recently developed FEM transfer matrix (or lead field bases) approaches [Weinstein et al., 2000; Gençer and Acar, 2004; Wolters et al., 2004] and advances in efficient FEM solver techniques for source analysis [Wolters et al., 2004] drastically reduce the complexity of the computations so that the main disadvantage of FEM modeling no longer exists. [Lanfer, 2007] compared run-time and numerical accuracy of a FEM source analysis approach (the FEM is based on a Galerkin approach applied to the weak formulation of the differential equation) using the Venant dipole model [Buchner et al., 1997] and the fast FE transfer matrix approach [Wolters et al., 2004] with a BE approach of [Zanow, 1997] (a double layer vertex collocation BE method [de Munck, 1992] using the isolated skull approach [Hämäläinen and Sarvas, 1989] and linear basis functions with analytically integrated elements [de Munck, 1992]) in combination with BE transfer matrices in an isotropic three layer sphere model. The reported numerical errors of the FE approach for realistic eccentricities in an isotropic three compartment sphere model were in the same range than those of the BEM approach, while, at the same time, the FE forward computation was faster than the BE forward computation. Additionally, similar errors and run-times were achieved with the FE approach in anisotropic four compartment sphere models, showing the large flexibility of this approach.
In this paper, we propose a low resolution conductivity estimation (LRCE) method using simulated annealing optimization in a realistically-shaped four compartment (scalp, skull, CSF and brain) finite element volume conductor model that individually optimizes the brain and the skull conductivity parameters while keeping the CSF conductivity fixed to the measurement of [Baumann et al., 1997] and the scalp conductivity to the value that is commonly used in source analysis [Buchner et al., 1997; Fuchs et al., 1998; Zanow, 1997; Waberski et al., 1998; Huiskamp et al., 1999]. The LRCE method uses a geometric model, in this case based on T1-/PD-MRI, and evoked potential data with high signal-to-noise ratio (SNR) as input. The method then determines the best combination of sources within a predefined source space together with the two individually optimized brain and skull conductivity values over a discrete parameter space, i.e., for each source and for each tissue conductivity the user has to define a reasonable set of a priori values. We evaluate the accuracy of the LRCE method in simulation studies before applying it to tactile somatosensory evoked potentials (SEP) with the focus on establishing the best values for the individual brain and skull conductivity. Besides using our new method for an improved source analysis of, e.g., SEP or auditory evoked potentials (AEP) (i.e., EEG data with a rather simple underlying source structure and a well-controlled and high SNR), the major future perspective for the LRCE is to provide an individually optimized volume conductor model (by means of exploiting the SEP/AEP data) that can then be used in a second step for the analysis of clinical or cognitive EEG data of the same subject/patient.
In the considered low frequency band (frequencies below 1000 Hz), the capacitive component of tissue impedance, the inductive effect and the electromagnetic propagation effect can be neglected so that the relationship between bioelectric fields and the underlying current sources in the brain can be represented by the quasi-static Maxwell equation
with homogeneous Neumann boundary conditions at the head surface
and a reference electrode with given potential, i.e., (ref) = 0, where σ is the conductivity distribution, is the scalar electric potential, 0 is the primary (impressed) current, Ω the head domain, Γ its surface and the surface normal at Γ [Plonsey and Heppner, 1967; Sarvas, 1987]. The primary current is generally modeled by a mathematical dipole at position 0 with the moment 0, 0 = 0 δ(−0) [Sarvas, 1987]. For a given primary current and conductivity distribution, the potential can be uniquely determined for what is known as the bioelectric forward problem.
For the numerical approximation of equations (1) and (2) in combination with the reference electrode, we use the finite element (FE) method. Different FE approaches for modeling the source singularity are known from the literature: a subtraction approach [Bertrand et al., 1991], a partial integration direct method [Weinstein et al., 2000], and a Venant direct method [Buchner et al., 1997]. In this study we used the Venant approach based on comparison of the performance of all three in multilayer sphere models, which suggested that for sufficiently regular meshes, it yields suitable accuracy over all realistic source locations [Wolters et al., 2007a, 2007b; Lanfer, 2007]. This approach has the additional advantage of high computational efficiency when used in combination with the FE transfer matrix approach [Wolters et al., 2004]. We used standard piecewise linear basis functions ϕi() =1 for = i, where i is the i-th FE node, and ϕj() = 0 for all j ≠ i. The potential is projected into the FE space, i.e., , where N is the number of FE nodes. Standard variational and FE techniques for equations (1) and (2) yield the linear system
where A is the stiffness matrix with dimension N × N,u the coefficient vector for h(N ×1), JVen the Venant approach right-hand-side vector (N ×1) [Buchner et al., 1997; Wolters et al., 2007a], and <·,·> the scalar product. A key feature of this study was to pursue solutions that achieve high computational efficiency. Let us assume that the S electrodes directly correspond to FE nodes at the surface of the head model (otherwise, interpolation is needed). It is then easy to determine a restriction matrix B (S−1)×N, which has only one non-zero entry with the value 1 in each row and which maps the potential vector u onto the potential vector at the (S−1) non-reference EEG electrodes, . With the following definition of the (S−1)×N transfer matrix for the EEG, T:= B A−1, a direct mapping of an FE right-hand side vector JVen onto the unknown electrode potentials is given. It was shown in [Wolters et al., 2004] how the transfer matrix T can efficiently be computed using an algebraic multigrid preconditioned conjugate gradient (AMG-CG) method. Note that JVen has only C non-zero entries (with C being the number of neighbors of the closest FE node to the source) so that TJVen only amounts in 2·(S−1)·C operations. Thus the resulting combination of the transfer matrix approach with the Venant method leads to implementations that are especially efficient [Lanfer, 2007], an essential feature for our study as will become clear in Section 2.3.
The non-uniqueness of the EEG inverse problem requires a combination of a viable forward problem, anatomical information, and a priori constraints on some aspect(s) of the solution. Here, we followed a dipole fit procedure that restricted the number of active sources to an application dependent number, K, of some few dipoles [Scherg and von Cramon, 1985; Mosher et al., 1992]. In addition, we defined an influence space with R discrete permissable source locations that was constrained to lay within the cortical gray matter. Given this influence space, the S scalp electrode locations, and a fixed volume conductor, we used the fast FE forward computation methods from Section 2.1 to compute a lead field matrix, L, which mapped sources directly to electrode potentials:
where J is a current source vector of dimension 3R×1 because we do not use the normal constraint, i.e., sources at the discrete source space nodes can have orientations in any direction. Φsim is the simulated potential vector of dimension S ×1 and L has dimension S×3R.
Since the potential depends linearly on the source moment (dipole orientation and strength) and nonlinearly on the source location, we use a two phase approach for source analysis [Buchner et al., 1997; Wolters et al., 1999]. We start with K initial source locations that are proposed by the non-linear optimization procedure simulated annealing (SA, see Section 2.2.2) and apply a linear least squares fit to the EEG data that determines uniquely the linear source orientation and strengths parameters, Jr(3K×1). The numerical solver for the linear least squares procedure employed a truncated singular value decomposition for the minimization [Wolters et al., 1999], based on a cost function, gf, that is the L2 norm of the difference between the simulated potential, Φsim(S×1), and the measured EEG potential, ΦEEG(S×1):
In this equation, Lr(S×3K) indicates the reduced lead field matrix for the current choice of source locations r = (1,···, K) with k the k-th source location (1 ≤ k ≤ K).
Since the volume conduction properties are incorporated in the lead field matrix Lr, the free nonlinear optimization parameters in this case are only the source locations. There is the choice between local optimization methods such as the Nelder-Mead simplex approach [Nelder and Mead, 1965] or the Levenberg-Marquardt algorithm [Marquardt, 1963] and global optimization approaches such as simulated annealing (SA) from combinatorial optimization [Kirkpatrick et al., 1983] or genetic algorithms [Kjellström, 1996]. In our paper, we decided for SA optimization because the challenge of local optimizers lies in determining the initial estimation of multiple parameters in the presence of local minima, and the global SA optimizer, often used when the search space is discrete as in our study, is generally more effective in localizing multiple parameters because it eliminates the need for high quality initial estimates [Haneishi et al., 1994; Gerson et al., 1994; Buchner et al., 1997; Uutela et al., 1998; Wolters et al., 1999]. A stochastic Metropolis acceptance test prevents the SA search from getting trapped in local minima as long as the cooling schedule is slow enough [Metropolis et al., 1953; Geman and Geman, 1984; Hütten, 1993]. For the cooling schedule, a so-called temperature t regulates the acceptance probability. Throughout the optimization process, t decreases according to a cooling rate ft. Initially, t is set to a high value, resulting in the acceptance of most new parameters (with even larger gf) and as the temperature decreases by means of ft, it is less likely for new parameters (with larger gf) to be accepted. This enables the search to focus on the vicinity of the minima at the later stages of the optimization process.
The proposed LRCE method adds electrical tissue conductivities as additional optimization parameters to the cost function to the already parameterized source locations. Here the set of optimization parameters including the conductivities was
where L is the number of tissue compartments and σl is the conductivity parameter for the l-th tissue compartment (1 ≤ l ≤ L). Each source location k was allowed to vary within the defined discrete influence space as described in Section 2.2. The conductivity σl of tissue compartment l was allowed to have its value from a predefined discrete set of possible conductivity values
Here, Hl is the number of possible conductivity values for tissue compartment l. We could choose Hl to be a large number (high resolution) for tissue l, but this would strongly increase computational costs and might be rather unrealistic given the limited SNR in measured EEG data. Therefore, we confined it to a rather small set of conductivity values (e.g., the different measured and estimated values for the considered head tissue that can be found in the literature).
Given the influence source space and the electrode locations, we precomputed a set of lead field matrices and collected them in Λ, which corresponded to all possible combinations of conductivity values for all tissue compartments of interest. This resulted in the number of lead field matrices in Λ.
with L(σh1, ···, σhL) being the (S× 3R) lead field matrix for the specific choice of conductivities. During each iteration of the SA method, the set of optimization parameters includes not just a new estimate of the bioelectric source, but a new configuration of both sources and conductivities in which we allow changing the value of only one parameter chosen randomly per iteration. By limiting the choice of conductivities to a discrete set of values, we maintain computational efficiency by applying the associated precomputed lead field matrix from the set Λ. The total number of possible configurations for sources and conductivities is
The SA optimizer searches for an optimal configuration of dipole source locations r = (1, ···, K) and tissue conductivities σ = (σ1, ···, σL) that ensure the best fit to the measured data:
The following summarizes the general procedure of the LRCE:
To carry out the LRCE analysis requires the construction of detailed realistic head models, in this case from MR image data. Here we outline the steps for constructing such a model. Our approach emphasizes accurate modeling of the CSF and skull compartments [Cuffin, 1996; Huiskamp et al., 1999; Ramon et al., 2004; Wolters et al., 2006; Wendel et al., 2008; Rullmann et al., 2008]. The influence of the skull thickness is closely related to the influence of skull conductivity and therefore especially important for a successful application of the presented LRCE algorithm [Cuffin, 1996; Huiskamp et al., 1999; Ramon et al., 2004]. To achieve the required accuracy of the head models, we made use of a combination of two different MRI modalities applied to a single subject. T1-weighted MRI is well suited for the segmentation of tissue boundaries like gray matter, outer skull and scalp. In contrast, the identification of the inner skull surface (defining thicknesses of skull and CSF compartment) is more successful from a Proton density MRI (PD-MRI) sequence because the difference in the quantity of water protons between intra-cranial and bone tissues is large. T1- and PD-MR imaging of a healthy 32 year-old male subject was performed, the images were aligned and segmented in a realistic four compartment (scalp, skull, CSF, brain) volume conductor model with special attention to the poorly conducting human skull and the highly conductive CSF following the procedures described in [Wolters et al., 2006]. The T1 images provided the information on soft tissues while the registered PD image enabled the segmentation of the inner skull surface and thus a correct modeling of skull and CSF compartmental thickness. In source reconstruction, it is generally accepted that the weak volume currents outside the skull and far away from the EEG sensors have a negligible influence on the measured fields [Buchner et al., 1997; Fuchs et al., 1998]. We therefore did not make any effort to segment the face and used instead a cutting procedure typical in source analysis based on realistically-shaped volume conductor modeling [Buchner et al., 1997; Fuchs et al., 1998].
Figure 1 shows the results of this approach for the segmentation of the inner skull/outer CSF surface compared with results from an estimation procedure that used exclusively the T1-MRI. The estimation procedure started from a segmented brain surface and estimated the inner skull by means of closing and inflating the brain surface.
A prerequisite for FE modeling is the generation of a mesh that represents the geometric and electric properties of the head volume conductor. To generate the mesh, we used the CURRY software [CURRY, 2000] to create a surface-based tetrahedral tessellation of the four segmented compartments. The procedure exploited the Delaunay-criterion, enabling the generation of compact and regular tetrahedra [Buchner et al., 1997; Wagner et al., 2000] and resulted in a finite element model with N= 245,257 nodes and 1,503,357 tetrahedra elements.
The FE mesh is shown in Figure 2.
An influence source space that represented the brain gray matter in which dipolar source activities occur was extracted from a surface 2 mm beneath the outer cortical boundary. The influence space was tessellated with a 2 mm mesh resulting in R = 21,383 influence nodes (shown in Figure 2). Since the influence mesh is only a rough approximation of the real folded surface and does not appropriately model the cortical convolutions and deep sulci, no normal-constraint was used, i.e., the dipole sources were not restricted to be oriented perpendicular to the source space. Instead, dipole sources in the three Cartesian directions were allowed.
Simulation studies were carried out to validate the new LRCE approach. For the reference FE volume conductor, isotropic conductivity values of 0.33 (see [Haueisen, 1996] and references therein), 0.0132 [Lai et al., 2005], 1.79 [Baumann et al., 1997], and 0.33 S/m (see [Haueisen, 1996] and references therein) were assigned to the scalp, skull, CSF, and brain compartment of the FE model from Section 3.2, respectively. This led to a brain:skull conductivity ratio (four-compartment head model) of 25 for the reference volume conductor. For the modeling of the EEG, 71 electrodes were placed on the reference volume conductor surface according to the international 10/10 EEG system. Two reference dipole sources were positioned on influence nodes in area 3b of the primary somatosensory cortex (SI) in both hemispheres, as shown in Figure 2 (right). Two source orientation scenarios were considered, in which both sources were either oriented quasi-tangentially or quasi-radially with regard to the inner skull surface. In both scenarios, the two sources were simultaneously activated using current densities of 100 nAm. Another experiment consisted of just a single source in the left SI with quasi-tangential or quasi-radial direction and a source strength of 100 nAm. Forward potential computations were carried out for the different scenarios using the direct FE approach as described in Section 2.1. Noncorrelated Gaussian noise was then added with SNR’s of 40, 25, 20, and 15 dB (SNR(dB):= 20 log10(SNR) with , where is the noisy signal and ε[i] the noise at electrode i).
Figure 3 shows the potential maps for the two-sources experiment for both orientation scenarios, the quasi-tangential (top row) and the quasi-radial orientations (bottom row) for different SNR values.
For the SA optimization, the source space from Section 3.2 was used as the influence space. A very slow cooling schedule with the cooling rate ft of 0.99 was applied in order to make sure that the search reached the global minimum of the cost function. The localization error was defined as the Euclidian distance between the somatosensory reference source locations and the inversely fitted ones resulting from the LRCE. The residual variance ν of the goal function was calculated as the percentile misfit between the noisy reference potential and the fitted potential that was computed from the fitted source parameters and conductivities. The explained variance shown in the result tables is 100% −ν.
We measured somatosensory evoked potential (SEP) data in order to apply our LRCE approach to real empirical EEG data. Tactile somatosensory stimuli were presented to the right index finger of the subject from Section 3.1 using a balloon diaphragm driven by bursts of compressed air. We compensated for the delay between the electrical trigger and the arrival of the pressure pulse at the balloon diaphragm as well as the delay caused by the inertia of the pneumatic stimulation device (half-way displacement of the membrane), together 52 ms in our measurements. Following standard practice [Mertens and Lütkenhöner, 2000], the stimuli were presented at 1 Hz (±10% variation to avoid habituation effects). A 63 channel EEG (10% system) recorded the raw time signals for the SEP study. Two electrooculography (EOG) electrodes were furthermore used for horizontal and vertical eye movement control. The collection protocol consisted of three runs of 10 minutes each EEG data with a sampling rate of 1200 samples/sec using a real time low pass filter of 0–300 Hz. The BESA software [BESA, 2007] was then used for a rejection of noise-contaminated epochs (e.g., epochs containing eye movements detected by means of the EOG channels) and for averaging the non-contaminated epochs within each run (83% of 601 epochs for run 1, 90% of 605 epochs for run 2 and 89% of 602 epochs for run 3). In order to optimize the SNR, the SEP data were furthermore averaged over the 1579 non-contaminated epochs of the three runs. The data was measured with FCz as reference electrode. The baseline-corrected (from −35 ms to 0 ms pre-stimulus) averaged EEG dataset was filtered using a 4th order butterfly digital filter with a bandwidth of 0.1 to 45 Hz. When using the prestimulus interval between −20 ms and 0 ms for the determination of the noise level and the peak of the first tactile component at 35.3ms as the signal, we achieved a SNR of 24dB. Finally, by means of a channel-selection procedure (exclusion of 20 ipsilateral electrodes with poor SNR), we were able to even increase the SNR to 26.4 dB.
A butterfly- and a position-plot of the SEP data is shown in Figure 4.
All simulations and evaluations ran on a Linux-PC with an Intel Pentium 4 processor (3.2GHz) using the SimBio software environment (SimBio, 2008).
We performed the LRCE procedure as described in Section 2.3 with an inverse two-dipole fit on the discrete influence space, while additionally allowing skull and brain conductivity to vary as free discrete optimization parameters. The permitted brain conductivities (σbrain) were 0.12, 0.33 [Haueisen, 1996], and 0.48 S/m. For each brain conductivity, the skull conductivity (σskull) was allowed to vary so as to achieve brain:skull ratios (four-compartment head model) of 80, 40, 25, 15, 10, 8, and 5. The CSF conductivity remained fixed at 1.79 S/m [Baumann et al., 1997] and the scalp conductivity at 0.33 S/m [Haueisen et al., 1996; Fuchs et al., 1998; Huiskamp et al., 1999]. Because of the fixed conductivities, possible problems are avoided that are due to the ambiguity between source strength and overall conductivity. This resulted in a total of 21 conductivity configurations.
Following equation (4), the total number of possible source and conductivity configurations in this simulation was thus approximately 4.8 billion.
Table 1 contains the LRCE source localization and conductivity estimation results for the simulated reference EEG data and Table 2 the LRCE reconstruction errors in the corresponding dipole moments. As the tables show, besides appropriately reconstructing both sources, the LRCE was able to accurately select the reference conductivity values of the brain and the skull compartment in the cases of no noise (max.errors: 0mm loc., 0 degree orientation, 0% magnitude) and low noise (40 dB, max.errors: 3mm loc., 6 degree orientation, 18% magnitude). However, for the noisier data with an SNR of 25 or lower, neither the somatosensory sources nor the brain and the skull conductivity values could be reconstructed correctly. The wall clock time for setting up the global leadfield matrix Λ was 199 minutes. When averaging over all noise configurations and source orientation scenarios, the SA needed about 17 hours of computation time for finding the global optimum (results indicated in Tables 1 and and2).2). Much of it was access-time to the global leadfield matrix within the LRCE procedure.
In the second simulation, we first generated noise-free and noisy reference data for a single dipole source in the left somatosensory cortex and then performed a single dipole fit with skull and brain conductivity as two additional free optimization parameters in the LRCE. We used the same scalp, skull, CSF, and brain conductivity values as in the previous simulation:
The number of possible source and conductivity configurations was 449K (equation (4)).
As shown in Tables 3 and and4,4, the conductivities were accurately estimated for reference data with 40dB and 25dB SNR and the source reconstruction errors were very low (max.errors for 40dB: 0mm loc., 1 degree orientation, 1% magnitude; max.errors for 25dB: 2mm loc., 4 degree orientation, 2% magnitude). For 20dB, the skull to brain conductivity ratio was still correct and the source reconstruction was still acceptable (max.errors: 4mm loc., 9 degree orientation, 12% magnitude), but the brain conductivity was no longer correctly reconstructed. Still higher noise levels led to unacceptable results. Like in Section 4.1.1, the wall clock time for setting up the global leadfield matrix Λ was 199 minutes. When averaging over all noise configurations and source orientation scenarios, the LRCE procedure took about 1.3 minutes of computation time for finding the global optimum (results indicated in Tables 3 and and44).
We carried out a third simulation, in which only skull conductivity was allowed to vary with fixed conductivity values for brain (0.33 S/m), scalp (0.33 S/m), and CSF (1.79 S/m). The brain:skull conductivity ratio (four-compartment model) was chosen as follows.
The total number of possible source and conductivity configurations for this scenario was 1.6 billion (equation (4)).
As shown in Tables 5 and and6,6, for both source orientation scenarios, the LRCE estimated the skull conductivity correctly down to a 20 dB SNR, while reasonable source reconstructions were only achieved down to 25 dB (<8mm loc., <16degree orientation, <15% magnitude). The LRCE reconstruction failed to give acceptable results for both the brain:skull conductivity ratio (four-compartment model) and the source reconstructions only at an SNR of 15dB or lower. The wall clock time for setting up the global leadfield matrix Λ was 66 minutes. When averaging over all noise configurations and source orientation scenarios, the LRCE procedure took about 451 minutes of computation time for finding the global optimum (results indicated in Tables 5 and and6).6). Again, much of it was access-time to the global leadfield matrix within the LRCE procedure.
In a last simulation, volume conductors with fixed skull conductivity values from the set of σskull from Section 4.1.3 were used. For these fixed volume conductors, only the two somatosensory sources were reconstructed on the discrete influence space using the simulated annealing optimizer with reference EEG data at an SNR of 25dB.
The results in Table 7 show the effects of an erroneous choice of the brain:skull conductivity ratio (four-compartment model) (80, 40, 15, 10, 8, 5) on the localization accuracy in comparison to the localization errors caused just by the addition of noise when using the correct brain:skull ratio (four-compartment model) of 1:25. Incorrect skull conductivity within the source localization caused large localization errors. As expected, the correct skull conductivity (σbrain/σskull = 25) gave the smallest localization errors and the highest explained variance for both source orientation scenarios.
In a last examination, the new LRCE algorithm was applied to the post stimulus P35 component of the averaged SEP data at the peak latency of 35.3ms as indicated in Figure 4. The detailed four compartment (scalp, skull, CSF, and brain) finite element model with improved segmentation of skull and CSF geometry described in Section 3.2 was used as the volume conductor. Because of the limiting SNR of 26.4 dB for the SEP data and based on our simulation results from Section 4.1, we focused on the simultaneous reconstruction of the contralateral somatosensory P35 source in combination with the estimation of both the brain and the skull conductivities. Accordingly, we assigned fixed isotropic conductivities to CSF (1.79 S/m) [Baumann et al., 1997] and scalp (0.33 S/m) [Haueisen et al., 1996; Fuchs et al., 1998; Huiskamp et al., 1999]. Again, the source space from Section 3.2 was used as the influence space for simulated annealing optimization together with brain:skull conductivity ratios (four-compartment head model) of 140, 120, 100, 80, 72, 60, 42, 25, 23, 15, 10, 8 and 5 ([Hoekema et al., 2003], who claimed ratios of 10 up to only 4).
The total number of possible source and conductivity configurations was 1,026K.
Applying the LRCE approach resulted in the contralateral somatosensory source shown in Figure 5, in the brain conductivity of 0.48S/m and in a skull conductivity of 0.004 S/m, with an explained variance of 99%. While the value of skull conductivity is close to what is generally used in three-compartment head model based source analysis, with 0.48S/m, the value of brain conductivity is higher than the commonly used (in three-compartment approaches) value of 0.33S/m (e.g., [de Munck and Peters, 1993; Buchner et al., 1997; Fuchs et al., 1998; Zanow, 1997; Waberski et al., 1998; Huiskamp et al., 1999]). The estimated brain conductivity is however still in the range of brain conductivity values that were determined by others (e.g., the value of 0.57S/m for subject 1 in [Goncalves et al., 2003a], 0.43S/m for subject 5 in [Goncalves et al., 2003b], 0.42S/m for subject S2 in [Baysal and Haueisen, 2004]). The wall clock time for setting up the global leadfield matrix Λ was about 315 minutes and the LRCE procedure took about 10 minutes of computation time.
We developed a low resolution conductivity estimation (LRCE) procedure to individually optimize a volume conductor model from a human head with regard to both geometry and tissue conductivities. We only exploited somatosensory evoked potential (SEP) data and a combined T1-/PD-MRI dataset for the construction of a four-tissue (scalp, skull, cerebrospinal fluid (CSF), brain) volume conductor FE model. The proposed procedure is safe and noninvasive, and EEG laboratories should most often have access to such datasets, so that no additional hard- and software is needed, in contrast to, e.g., approaches based on Electrical Impedance Tomography (EIT) [Gonçalves et al., 2003b]. For the FE model, a special focus was on an improved modeling of the skull shape and thickness and on the highly conducting CSF compartment [Baumann et al., 1997; Huiskamp et al., 1999; Ramon et al., 2004; Wolters et al., 2006; Wendel et al., 2008; Rullmann et al., 2008]. Obtaining accurate skull geometry is important because changes in skull conductivity are known to be closely related to changes in its compartmental thickness. The correction for geometry errors in modeling the skull compartment were furthermore shown to be essential for the measurement of skull conductivity [Gonçalves et al., 2003b]. While other authors have used parameter estimation in continuous parameter space with local optimization algorithms [Fuchs et al., 1998; Gutiérrez et al., 2004; Vallaghe et al., 2007, Zhang et al., 2006], we propose the combination of a discrete low resolution parameter estimation with a global optimization method applied to realistic four-compartment geometry to better take into account the limited signal-to-noise (SNR) of real SEP or auditory evoked potential (AEP) measurement data. Because the cost function is shallow [Gonçalves et al., 2003a], the proposed procedure using realistic FE volume conductor modeling and simulated annealing (SA) optimization for approximating the global minimum in acceptable computation time is important. While other authors used three compartment boundary element (BE) [Fuchs et al., 1998; Gonçalves et al., 2003a;Plis et al., 2007; Vallaghe et al., 2007] or finite element models [Zhang et al., 2006] (in the latter, additionally to the three layers scalp, skull and brain, a low conducting silastic ECoG grid was modeled) for conductivity estimation, we additionally model the CSF with a fixed conductivity of 1.79S/m [Baumann et al., 1997], not only because its modeling was shown to have a large impact on forward and inverse source analysis [Huang et al., 1990; Baumann et al., 1997; Ramon et al., 2004; Wolters et al., 2006; Wendel et al., 2008; Rullmann et al., 2008], but also to avoid the problem of the ambiguity between source strength and overall conductivity. In [Rullmann et al., 2008], non-invasive EEG source analysis was validated by means of intra-cranial EEG measurements and it was shown that ignoring the CSF by means of the commonly used three-compartment realistically-shaped volume conductor led to spurious reconstruction results. [Plis et al., 2007] derived a lower Cramer-Rao bound for the simultaneous estimation of source and skull conductivity parameters in a sphere model for dipoles whose locations were not constraint within the inner sphere volume. Since source depth and skull conductivity are closely related, their final result was that it is impossible to simultaneously reconstruct both source and skull conductivity parameters from measured surface EEG data in the sphere model. This is an important theoretical result, however, there are strong differences to our study. Our study, as well as the symmetric BEM study of [Vallaghe et al., 2007], used a cortex constraint, i.e., sources were only allowed on a surface. We furthermore used a realistic four-compartment FE model of the head instead of the spherical volume conductor model that was used for the derivation of the Cramer-Rao bounds in [Plis et al., 2007] and we fixed the conductivity of the CSF compartment in our analysis to the value measured by [Baumann et al., 1997]. We only allowed a user-given discrete set of some few (“low resolution”) possible conductivity values for those tissues where conductivity measurements or other methods resulted in different estimates. We propose to only apply the presented LRCE algorithm to EEG data where the underlying sources are rather simple and where very good SNR ratios can be achieved like, e.g., SEP and/or AEP data.
In the first simulation studies, we evaluated the LRCE algorithm in EEG simulations for its ability to determine both the brain and the skull tissue conductivities together with the reconstruction of one and two reference sources. At relatively low noise levels (down to 25 dB SNR in the single source scenario and down to 40 dB SNR in the two source scenario), the LRCE resulted in acceptable reconstruction errors for the reference sources and correctly estimated reference tissue conductivities, while results became unstable when further increasing the noise. We also set up a simulation for the reconstruction of the skull to brain conductivity ratio (four-compartment model) together with two sources in which results were reasonable (correct skull:brain conductivity ratio, max. source reconstruction errors: <8mm loc., < 16degree orientation, <15% magnitude) down to noise levels of 25 dB. We found in our simulations that the most accurate source reconstructions were associated with the correctly estimated conductivities (or conductivity ratio) and, moreover, that assuming an incorrect conductivity ratio had a profoundly negative effect on the source reconstruction accuracy.
In a last examination, we applied the LRCE to measured tactile SEP data with the focus on estimating both the brain and the skull conductivity. With an SNR of 26.4 dB, the data were in the noise range of the second simulation study, which was based on a single equivalent current dipole model. As shown in numerous studies [Mertens and Lütkenhöner, 2000; Hari and Forss, 1999], this source model is adequate because the early SEP component arises from area 3b of the primary somatosensory cortex (SI) contralateral to the side of stimulation. Our explained variance to the measured data of about 99% for this source model further supports our choice. The results from the LRCE analysis were a brain conductivity of 0.48 S/m and a skull conductivity of 0.004 S/m. While this skull conductivity corresponds to the traditional value in the literature [de Munck and Peters, 1993; Buchner et al., 1997; Fuchs et al., 1998], we found the brain to have a lower resistance than generally assumed in three-compartment head modeling approaches (e.g., [de Munck and Peters, 1993; Buchner et al., 1997; Fuchs et al., 1998; Zanow, 1997; Waberski et al., 1998; Huiskamp et al., 1999]), but it is however still in the range of brain conductivity values that were determined by others (e.g., the value of 0.57S/m for subject 1 in [Goncalves et al., 2003a], 0.43S/m for subject 5 in [Goncalves et al., 2003b], 0.42S/m for subject S2 in [Baysal and Haueisen, 2004]).. Many recent papers have focused on the brain:skull conductivity ratio and a large variability of results have been reported for this value including 80 [Homma et al., 1995], 72 [Gonçalves et al., 2003a], 42 [Gonçalves et al., 2003b], 25 ± 7 [Lai et al., 2005], 23 [Baysal and Haueisen, 2004], 18.7 ± 2.1 [Zhang et al., 2006], 15 [Oostendorp et al., 2000] and 8 [Hoekema et al., 2003]. Because of the higher conductivity of the brain, with an estimated brain to skull conductivity ratio of 120 (in a four-compartment head model), our LRCE result is larger than the commonly used ratio of 80 [Rush and Driscoll, 1968,1969; Homma et al., 1995]. Note that, in contrast to the studies using three-compartment modeling [Homma et al., 1995; Oostendorp et al., 2000; Gonçalves et al., 2003a; Gonçalves et al., 2003b; Baysal and Haueisen, 2004; Lai et al., 2005; Vallaghe et al., 2007], our approach took the highly conducting CSF compartment into account. It is well known that an increased conductivity of the brain compartment leads to a decreased potential magnitude at the head surface while an increased conductivity of the CSF leads to an increased potential magnitude. Since we modeled the CSF with the value of 1.79S/m as measured by [Baumann et al., 1997] (i.e., a more than a factor of 5.4 higher value than the commonly used 0.33S/m in the three-compartment models), an increased conductivity value for the brain compartment has to be expected in the four-compartment model. Our brain-to-skull conductivity ratio (in a four-compartment model) of 120 thus has to be interpreted in light of the above considerations.
With regard to computational complexity (or feasibility in daily routine), former FE approaches, which were not based on the presented transfer matrix approach and on algebraic multigrid FE solver methods would have needed weeks or even months for the computation of a single leadfield matrix for a single conductivity configuration so that the proposed FE-based LRCE approach would not have been feasible in practice. In [Buchner et al., 1997], the computation of a single leadfield matrix for an FE mesh with 18,322 nodes and an influence space with 2,914 nodes took roughly a week of computation time. [Waberski et al., 1998] used an FE model with 10,713 nodes and concluded that improved headmodeling by finer discretization and more accurate representation of the conductivities are necessary and parallel computing is needed to speed up the computation. The FE head model of [Zhang et al., 2006] for the estimation of the in-vivo brain-to-skull conductivity ratio had 29,858 nodes. For our presented LRCE approach, the underlying FE mesh had a resolution of 245,257 FE nodes, which was necessary not only to appropriately model the CSF compartment, and our influence space had 21,383 nodes. Furthermore, our LRCE algorithm does not only need to precompute a single leadfield matrix, but as many leadfield matrices as we have combinations of user-given conductivity values for the different tissue compartments as indicated by the global leadfield matrix in equation (3). Still the presented LRCE approach, as indicated by means of the computation times in Section 4 (measured on a single processor machine, see Section 3.5), is practically feasible in daily routine with a computational amount of work in the range of some few hours.
The following limitations of our study are important: The data of a single subject is not representative for other subjects since we have to be aware of larger inter- and intra-subject variability. The variability can be related to age, diseases, environmental factors, and personal constitution as shown in animal studies [Crile et al., 1922] and as shown for humans by means of the large discrepancy in the estimated brain-to-skull conductivity ratios (in three-compartment models) between 80 [Rush and Driscoll, 1968, 1969; Homma et al., 1995] and 15 [Oostendorp et al., 2000]. Further simulation studies should be carried out that consider noise from, e.g., the pre-stimulus interval of evoked potential measurements. The presented LRCE procedure has to be automatized in order to allow a statistical evaluation of possible errors and instabilities at different noise levels. We are currently working on such investigations for a combined SEP/SEF-LRCE approach. The influence of a realistic extent of an active cortical patch on our focal-source based LRCE method should be evaluated and its sensitivity to biological noise (non-modeled “noise” current sources in the brain) has to be examined. The performance of other global optimization approaches such as genetic algorithms [Kjellström, 1996] should be compared with the approach chosen here and higher FE resolutions have to be used in order to avoid geometry representation problems in areas where, e.g., the CSF or the skull compartments are very thin using, e.g., 1mm hexahedra FE modeling as described in [Rullmann et al., 2008].
The current results illustrate the feasibility of building an optimized volume conductor model with regard to both geometry and conductivity. As we have formulated it, such a study requires accurate head geometry, in this case from both T1- and PD-weighted MRI (or T2-MRI) and cortical constraints on the sources. The highly conducting CSF should not be neglected in the headmodel [Huang et al., 1990; Baumann et al., 1997; Ramon et al., 2004; Wolters et al., 2006; Wendel et al., 2008; Rullmann et al., 2008] and our procedure takes this compartment into account. By obtaining SEP data, which allows independent reconstruction of the underlying bioelectric source, it is then possible to estimate the optimal conductivities for the individual subject using the proposed LRCE approach in highly realistic finite element models, provided that the data has a sufficient SNR ratio. Note that one might also think of simultaneously evaluating SEP data of different finger or toes ([Vallaghe et al., 2007], e.g., used left and right hand index finger SEP data). A related finding from this study is, there is a trade off between the number of independent parameters that can be determined and the complexity of the assumed source model. The specific trade off point is also strongly influenced by the quality of the measured electric potentials. Thus the number of parameters that can be dependably estimated is a function of both the signal quality and the number and quality of a priori knowledge about, for example, the source location or orientation through a combination with fMRI or anatomical and/or functional arguments (e.g., a strong restriction of the source location to anatomically and physiologically reasonable areas close to the somatosensory SI area). In this context, others have suggested that by including MEG data in the scheme [Fuchs et al., 1998; Huang et al., 2007], it will be possible to improve stability considerably. We note that our approach differs from their procedures with regard to both head modeling and conductivity optimization.
The success of the conductivity optimization approach and the more general advantages of customized geometric models suggest a procedure for clinical applications. First of all, one could use SEP and/or AEP data with high SNR together with T1- and PD-MR (or T2-MR) images from the patient to construct a model that would be optimized for both geometric accuracy and individual conductivity values. With this volume conductor model in place, recorded potentials from more complex and clinically interesting sources could drive the inverse solution and source analysis.
A better approximation to the real volume conductor using the proposed LRCE method is an important step towards simultaneous EEG/MEG source analysis. Combining EEG and MEG modalities compensates each others disadvantages, i.e., poor sensitivity of MEG to radial sources and the much stronger conductivity dependency of EEG [Fuchs et al., 1998; Huang et al., 2007]. Using combined somatosensory evoked potentials and fields (SEP/SEF) in combination with T1-and PD-MRI (or T2-MRI) should further stabilize the application of the presented LRCE method for the estimation of tissue conductivities. For the quasi-tangentially oriented P35 somatosensory source, MEG-SEF data can be exploited to strongly restrict the source location and especially its depth as shown, e.g., in [Fuchs et al., 1998; Huang et al., 2007], so that the resolution of the proposed LRCE method with regard to the conductivities of the different compartments could be increased. With such data in hand, the presented LRCE method using FE volume conductor modeling might also contribute to the estimation of conductivity values for further compartments like the scalp or of anisotropy ratios in the skull and brain compartments [Marin et al., 1998; Haueisen et al., 2002; Wolters et al., 2006; Rullmann et al., 2008].
This research was supported by the Center for Integrative Biomedical Computing, NIH NCRR Project 2-P41-RR12553-07 and by the German Research Foundation (DFG), projects WO1425/1-1 and JU 445/5-1. The authors would like to thank the anonymous reviewers for their helpful critics and comments that significantly improved our manuscript.