Cortical potential imaging using L-curve and GCV method to choose the regularisation parameter
© Subramaniyam et al; licensee BioMed Central Ltd. 2010
Published: 3 June 2010
The electroencephalography (EEG) is an attractive and a simple technique to measure the brain activity. It is attractive due its excellent temporal resolution and simple due to its non-invasiveness and sensor design. However, the spatial resolution of EEG is reduced due to the low conducting skull. In this paper, we compute the potential distribution over the closed surface covering the brain (cortex) from the EEG scalp potential. We compare two methods – L-curve and generalised cross validation (GCV) used to obtain the regularisation parameter and also investigate the feasibility in applying such techniques to N170 component of the visually evoked potential (VEP) data.
Using the image data set of the visible human man (VHM), a finite difference method (FDM) model of the head was constructed. The EEG dataset (256-channel) used was the N170 component of the VEP. A forward transfer matrix relating the cortical potential to the scalp potential was obtained. Using Tikhonov regularisation, the potential distribution over the cortex was obtained.
The cortical potential distribution for three subjects was solved using both L-curve and GCV method. A total of 18 cortical potential distributions were obtained (3 subjects with three stimuli each – fearful face, neutral face, control objects).
The GCV method is a more robust method compared to L-curve to find the optimal regularisation parameter. Cortical potential imaging is a reliable method to obtain the potential distribution over cortex for VEP data.
The electro- and magneto-encephalography (EEG, MEG) are the tools that are commonly used to record the brain activity on a sub millisecond time scale . One of the important goals of EEG/MEG research is the localisation of the underlying sources giving rise to the activities on the scalp . Such a localisation is not possible using the scalp EEG/MEG data alone. This stems from the theory of classical electromagnetic inverse problems, which says that different source configurations can give rise to same electric potential and magnetic field on the scalp .
The non-uniqueness in the inverse problem can be solved by determining the potential distribution over a closed surface, which encloses all the sources in the brain . Strictly speaking, this is not an inverse solution as no information about the generators is obtained. Also, the same potential distribution over the cortex can be explained by more than one configuration of generators. The advantages of computing cortical potential distribution are a) It can be determined uniquely as solution space is restricted to a closed surface  and b) It might provide better and enhanced picture of what is going on inside the brain with significant hints about the suppositional sources. This is a spatial enhancement technique and the scalp potential is de-blurred and computed over the cortex. There exist many spatial enhancement algorithms [5–7]. In cortical potential imaging technique a volume conductor model of the head is constructed (realistic or spherical). The properties of the region between the cortex and scalp are only considered, which has no sources  and the Laplace’s equation is solved.
The tissue conductivity is an important issue in volume conductor modeling. Over the years, different values for tissue conductivity have been reported for head model. A wide range of values for the conductivity of skull has been reported. In their seminal article, Rush and Driscoll reported that the resistivity of skull was 80 times that of brain . For a very long time, this value was used in modeling studies. In the year 2000 Ooestendorp et al.  reported a new value for resistivity ratio and it was 1:15:1 for brain:skull:scalp. Recent values for skull resistivity are reported in [11–13]. It is worth noting that, the accuracy of the inverse solution is affected by the value chosen for skull resistivity. Also, skull contains sinuses and areas of thinning. This impacts the volume current distribution and makes the current path longer in the areas which has holes. This makes localisation difficult .
The EEG inverse problem is ill-posed due to the fact that the transfer matrix, which relates the sources to measurements, is ill-conditioned. One way to deal with the ill-posedness is to apply regularisation techniques. There exist many regularisation techniques, the primary ones being Tikhonov regularisation and truncated singular value decomposition (TSVD). A key issue in the application of such regularisation methods is the selection of the optimal regularisation parameter. A detailed review about these methods can be found in .
Cortical potential distribution is of clinical importance and copious research has been done to enhance the spatial details of EEG [5–8]. By looking at the cortical potential distribution, a surgeon can do a better pre-surgical planning for e.g. in case of epileptic surgery . Such information is also valuable to psychologists in the field of cognitive neuroscience. The VEP distribution over the cortex, for instance can provide better hints than scalp potential about the parts of the cortex which are more active for a certain external visual stimuli (faces, objects etc) .
The purpose of this study was to apply cortical potential imaging technique to a realistic head model constructed using VHM dataset and compare the L-curve and GCV method used to find the regularisation parameter. We used a real set of EEG data to demonstrate the effect of choosing regularization parameters using these two methods on the potential distribution over cortex.
A realistic FDM model of the human head was constructed using VHM dataset, which excluded the brain and had cortex defined as closed surface. The same model was used for all the three subjects, without any modifications made taking the anatomical features of the subjects. The brain can be excluded based on Gauss’s law since the sources exist only inside the brain.
Forward transfer matrix
The FDM model consisted of 103799 cortical nodes. The cortical nodes were defined as the nodes lying at the interface of CSF and the empty space. In order to find a relation between cortical potential and the scalp potential, a forward transfer matrix was obtained. Since the number of cortical nodes is high, which would result in higher computational time, the nodes were grouped into equal sized source areas. In order to achieve this, first a sphere was divided into n equal sized source areas. This sphere was placed at the origin of the cortex and these areas were projected on the surface of cortex.
Superposition principle was used to obtain the forward transfer matrix. One source area at a time was made active (i.e. setting one source area to 1 volt), and the rest were set to zero . The resulting potential at electrode locations on scalp was found. This information forms one column of the transfer matrix. This procedure was repeated for the rest of the source areas, resulting in a matrix of size m x n, where m is the number of electrode locations and n, the number of source areas. The relationship between the number of source areas and number of reconstructable basis vectors for different noise levels has been investigated by Ryynänen et al  for spherical head model. Based on this study, the cortex was divided into 2000 equal sized source areas.
Actual EEG potentials
The data used in this study was obtained from the EEG databank from the Biomedical Engineering Laboratory of Tampere University of Technology (TUT), Finland. The measurements were carried out at the research laboratory of the Biomedical Engineering department of TUT at FinnMedi, Tampere, Finland. The EEG data used was the N170 component of VEP for three subjects S1, S2 and S3 under three different stimuli – fearful face, neutral face and control object. The three subjects had age range from 18 to 28 years and had no neurological history. It was a 256-channel data, with sampling rate of 1000 Hz and had 40 epochs. The EEGLAB  toolbox in MATLAB was used to average the potential data. The time period was from -100 milliseconds to 500 milliseconds. The signal at 170 ms was extracted from the averaged data for all the three subjects and was used to compute the inverse solution. The gain of the amplifier was set to 5000 and the signal was filtered with band pass filter of frequencies between 0.1 Hz to 200 Hz.
The realistic electrode locations obtained were projected on the scalp surface of the FDM model. The scalp surface nodes were extracted from the FDM model as the nodes lying at the interface of skin and empty space. The co-ordinates of the points obtained from the Polhemus system was transformed into a new orthogonal co-ordinate system defined using the image dataset. The nasion and the pre-auricular points were on the x-y plane. The nasion was along the x-axis and the pre-auricular points were on the y-axis. The realistic electrode locations were first projected on a unit sphere. This sphere was placed at the origin of the new orthogonal co-ordinate system and the electrode locations were projected back onto the scalp surface.
Discrete Picard condition
Determination of the decay rate of Fourier co-efficients with respect to the singular values is important to investigate the properties of regularised solution. The discrete Picard condition states that, if the Fourier co-efficients on an average, decay faster than the singular values of the transfer matrix, then the regularised solution will have the same regularity properties as the exact solution . To investigate this, we obtained the discrete Picard plot for the measurement data concerning all the three stimuli for subjects S1, S2 and S3. The regularisation toolbox was used to obtain Picard plots.
Estimation of the cortical potential distribution
The forward transfer matrix A, vector of measured scalp potential b and the solution vector x are related by,
Ax = b (1)
To obtain the inverse solution, zero-order Tikhonov regularisation was applied . It is given by,
X = (AT A + αI)–1AT b (2)
where the factor denotes the matrix that maps the known vector, b onto the regularized solution. The regularisation toolbox  was used to compute the parameter and apply regularisation technique.
Comparison of cortical potential distribution
The cortical potential distribution obtained using GCV and L-curve method can be compared by computing what is known as the MAG index and RDM index .
where, n is the number of source areas, XGCVi and XLi are the cortical potential values for i th source area computed using GCV and L-curve methods respectively.
The MAG index is used to compare the magnitude of potential values, where as the RDM index gives the fitting quality of spatial distribution of potential. In case of perfect fitting, the MAG would be 1 and RDM would be 0 .
The discrete Picard plot
The L-curve and GCV functional curve
Regularisation Parameter, MAG, RDM
Cortical potential distribution for N170
The selection of regularisation parameter is the key issue in the implementation of regularisation methods for solving ill posed problems. It is the regularisation parameter that controls the smoothness of the solution and tuning of this parameter is essential in arriving at a meaningful and reliable solution. Two widely used and popular methods, L-curve and GCV method are compared in this study.
The L-curve method has emerged as a popular method over the past few years . The L-curve method generally produces a sharp, convex corner shaped as ‘L’. However, this always may not be the case. The L-curve might breakdown due to the dominance of correlated geometry noise over the Gaussian measurement noise and such instances have been reported in literature  or due to the violation of the discrete Picard condition . The GCV method is based on successively leaving out elements of measurement vector and computing the regularised estimate of the solution vector from the reduced dataset. It predicts the left out response from this newly estimated model. Unlike L-curve, under the circumstances of violation of discrete Picard condition, GCV method has performed well [28, 29].
The EEG potential data acquired for subject S1 was not of good quality as the electrodes near the right ear had high impedances. It was observed that, the discrete Picard condition was not satisfied in case of S1 as seen from Figure 2, which may have lead to the poor performance and breakdown of L-curve. Such instances of breakdown have been reported in the past [27, 30].
In this study, we observed that in the L-curve for data concerning S1, the residual norm was monotonically increasing with no considerable change in the semi-norm of the solution as evident from Figure 3. The regularised solution obtained using the regularisation parameter from L-curve for S1 was over-smoothed and blurred. This can be seen in Figure 5. The GCV method performed relatively well in case of S1, in spite of the violation of discrete Picard condition. The GCV curve gave a sharp minimum for the fearful face stimuli. In case of neutral face and control object, the GCV method still performs well maintaining its characteristic shape and giving a well defined minimum. In case of S2 and S3, where the data was of good quality and the discrete Picard Condition was satisfied, both the methods gave comparable regularisation parameter, which in fact were not very different from each other.
We also computed the cortical potential distribution for the N170 component using the regularisation parameter obtained from this study for all the three subjects. There have been studies in the past regarding the source localisation of the N170 component to faces and objects [31–33], but none have focussed on the cortical potential distribution. These studies have demonstrated that the probable sources for facial perception are located mainly in the fusiform gyrus, posterior superior temporal sulcus region, predominantly on the right hemisphere. In this study, our aim was not to localise the sources, but to compute the equivalent of scalp potential distribution over a closed surface, defined as cortex. Though such a distribution gives no information about the underlying complex set of sources, but it can be helpful in providing hints about the sources. From the cortical potential distribution computed from this study for the N170 component of the ERP, we can see that the maximal negative potential occurs around the occipital and temporal parts of the cortex. The fusiform face area, which is thought to be responsible for facial recognition is located on the fusiform gyrus and is larger in the right hemisphere, though there is slight variation from person to person. The cortical potential distribution computed in this study seemed consistent with this fact by displaying areas of high negative potential on the right hemisphere of cortex and such a potential distribution probably may give significant hints about the putative sources. To our knowledge, there has been no study that has aimed to analyse the potential distribution over the cortex for the N170 component. Though many studies have aimed at source localization, computation of potential distribution over the cortex is a lesser explored research area.
In future work, we plan to use subject – specific models and automatically digitise the electrode locations. Also, an accurate representation of cortex is needed including all the sulcus and gyri, to get a better spatial detail for the potential distribution.
The GCV method is more robust method than the L-curve method in finding the regularisation parameter even when there is noise dominance and discrete Picard condition is not satisfied. Cortical potential imaging technique is a feasible method to investigate the cortical potential distribution for N170 component of the VEP data.
We would like to thank PhD Noriyuki Takano, for providing the finite different method software.
This article has been published as part of Nonlinear Biomedical Physics Volume 4 Supplement 1, 2010: Consciousness and its Measures: Joint Workshop for COST Actions Neuromath and Consciousness. The full contents of the supplement are available online at http://www.nonlinearbiomedphys.com/supplements/4/S1.
- Michel CM, Murray MM, Lantz G, Gonzalez S, Spinelli L: EEG source imaging. Clinical Neurophysiology. 2004, 115: 2195-2222. 10.1016/j.clinph.2004.06.001.View ArticleGoogle Scholar
- Ford RM, Sidman RD, Ramsey G: Spatio-temporal progression of AEP P300 component using the cortical imaging technique. Brain Topography. 1993, 6: 43-70. 10.1007/BF01234126.View ArticleGoogle Scholar
- Helmholtz HLF: Ueber einige Gesetze der Vertheilung elektrischer Stro¨mein ko¨rperlichen Leitern mit Anwendung aud die thierisch-elektrischen Versuche. Ann Physik und Chemie. 1853, 9: 211-33. 10.1002/andp.18531650603.View ArticleADSGoogle Scholar
- Yamashita Y: Theoretical studies on inverse problems of electrocardiography and the uniqueness of solution. IEEE Trans Biomed Eng. 1982, 9: 719-725. 10.1109/TBME.1982.325002.View ArticleGoogle Scholar
- Babiloni F, Babiloni C, Carducci F, Fattorini L, Onorati P, Urbano A: Spline Laplacian estimate of EEG potentials over a realistic magnetic resonance-constructed scalp surface model. Electroencephalogr ClinNeurophysiol. 1996, 98: 363-7. 10.1016/0013-4694(96)00284-2.View ArticleGoogle Scholar
- He B, Lian J, Li G: High-resolution EEG a new realistic geometry spline Laplacian estimation technique. Clin Neurophysiol. 2001, 112: 845-52. 10.1016/S1388-2457(00)00546-0.View ArticleGoogle Scholar
- Gevins AS, Le J, Brickett P, Reutter B, Desmond J: Seeing through the skull: advanced EEGs use MRIs to accurately measure cortical activity from the scalp. Brain Topogr. 1991, 4: 125-31. 10.1007/BF01132769.View ArticleGoogle Scholar
- Ollikainen JO, Vauhkonen M, Karjalainen PA, Kaipio JP: A new computational approach for cortical imaging. IEEE Trans on Medical Imaging. 2001, 20: 325-332. 10.1109/42.921481.View ArticleGoogle Scholar
- Rush S, Drischoll DA: EEG Electrode Sensitivity. An application of reciprocity. IEEE Transactions on Biomedical Engineering. 1969, 16: 15-22. 10.1109/TBME.1969.4502598.View ArticleGoogle Scholar
- Oostendorp TF, Delbeke J, Stegeman DF: The Conductivity of the human skull: Results of in vivo and in vitro measurements. IEEE Transactions on Biomedical Engineering. 2000, 36: 1487-1491. 10.1109/TBME.2000.880100.View ArticleGoogle Scholar
- Akhtari M, Bryant HC, Mamelak AN, Heller L, Shih JJ, Mandelkern M, Maltachov A, Ranken DM, Best ED, Sutherling WW: Conductivities of three layered human skull. Brain Topography. 2002, 14: 157-161. 10.1023/A:1014590923185.View ArticleGoogle Scholar
- Haueisen D, Tuch DS, Ramon C, Schimpfs PH, Wedeen VJ, George JS, Belliveau JW: The influence of brain tissue anisotropy on human EEG and MEG. Neuroimage. 2002, 15: 159-166. 10.1006/nimg.2001.0962.View ArticleGoogle Scholar
- Hoekema R, Wieneke GH, Leijten FSS, van Veelen CWM, Rijen PC, Huiskamp GJM, Ansems J, van Huffelen AC: Measurement of the conductivity of skull temporarily removed during epilepsy surgery. Brain Topography. 2003, 16: 29-38. 10.1023/A:1025606415858.View ArticleGoogle Scholar
- Regan D: Human brain electrophysiology – evoked potentials and evoked magnetic fields in science and medicine. Elsiever Science Publishing. 1989Google Scholar
- Grech R, Cassar T, Muscat J, Camilleri KP, Fabri SG, Zervakis M, Xanthopoulus P, Sakkalis V, Vanrumste B: Review on solving inverse problem in EEG source analysis. Journal of Neuroengineering and Rehabilitation. 2008, 5: 25-10.1186/1743-0003-5-25.View ArticleGoogle Scholar
- He B, Zhang X, Lian J, Sasaki H, Wu D, Towle VL: Boundary element method-based cortical potential imaging of somatosensory evoked potentials using subjects’ magnetic resonance images. NeuroImage. 2002, 16: 564-576. 10.1006/nimg.2002.1127.View ArticleGoogle Scholar
- Leppänen JM, Kauppinen P, Peltola MJ, Hietanen JK: Differential electrocortical responses to increasing intensities of fearful and happy emotional expressions. Brain Research. 2009, 1166: 103-109. 10.1016/j.brainres.2007.06.060.View ArticleGoogle Scholar
- Walker SJ, Klipatrick D: Forward and inverse calculations using resistor network models of human torso. Circulation Research-Journal of American Heart Association. 1987, 61: 504-513.Google Scholar
- Ryynänen O, Hyttinen J, Malmivuo J: Study on the spatial resolution of EEG – effect of electrode density and measurement noise. Proceedings of the 26th Annual International Conference of the IEEE EMBS: 1-5. 2004, September ; San FranciscoGoogle Scholar
- Delorme A, Makeig S: EEGLAB: An open source toolbox for analysis of single-trial EEG dynamics. Journal of Neuroscience Methods. 2004, 134: 9-21. 10.1016/j.jneumeth.2003.10.009.View ArticleGoogle Scholar
- Hansen PC: Advances in computational bioengineering. Computational inverse problems in electrocardiology. 2000, The L-curve and its use in the numerical treatment of inverse problemsGoogle Scholar
- Tikhonov AN: Solution of incorrectly formulated problems and the regularization method. Soviet Math Dokl. 1963, 1035-1038.Google Scholar
- Hansen PC: Rank deficient and discrete ill posed problems. Numerical aspects of linear inversion. SIAM, Philadelphia. 1997Google Scholar
- Golub GH, Van Loan CF: Matrix Computations. 2000, Baltimore John Hopkins University Press, 3rdGoogle Scholar
- Hansen PC: Regularization tools, A Matlab package for ananlysis and solution of discrete ill posed problems.Google Scholar
- Johnson PR, Guljarini R: A new method for regularisation parameter determination in inverse solution of electrocardiography. IEEE Transactions on Biomedical Engineering. 1997, 44: 19-39. 10.1109/10.553710.View ArticleGoogle Scholar
- Vogel CR: Non-convergence of the L-curve regularization parameter selection method. Inverse Problems. 1996, 12: 535-547. 10.1088/0266-5611/12/4/013.MATHMathSciNetView ArticleADSGoogle Scholar
- Wabha G: Practical approximate solutions to linear operator equations when the data are noisy. SIAM J. Numer. Annal. 1977, 14: 651-667. 10.1137/0714044.View ArticleADSGoogle Scholar
- Wabha G: Spline Models for Observational Datavijay. 1990, SIAMGoogle Scholar
- Rodriguez G, Theis D: An algorithm for estimating the optimal regularization parameter by the L-curve. Rendiconti di Mathematica. 2005, 25: 69-84.MATHMathSciNetGoogle Scholar
- Defke I, Sandar T, Heidenreich J, Sommer W, Curio G, Trahms L, Lueschow A: MEG/EEG sources of the 170 ms response to faces co-localized in the fusiform gyrus. NeuroImage. 2007, 35: 1495-1501. 10.1016/j.neuroimage.2007.01.034.View ArticleGoogle Scholar
- Itier JR, Taylor M: Source analysis of N170 to faces and objects. Brain Imaging. 2004, 15: 1261-1265.Google Scholar
- Kanwisher N, Moscovitch M: The cognitive neuroscience of face processing: An introduction. Cognitive Neuropsychology. 2007, 17: 1-11.View ArticleGoogle Scholar
- Chauveau N, Franceries X, Aubry F, Celsis P, Rigaud B: Cortical Imaging on a Head Template: A Simulation Study Using a Resistor Mesh Model (RMM). Brain Topography. 2008, 6: 52-60. 10.1007/s10548-008-0059-0.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.