Conjugate Gradient Method Applied to Cortical Imaging in EEG/ERP

X. Franceries1,2,4, N. Chauveau1,2,*, A. Sors3, M. Masquere4 and P. Celsis1,2 1Inserm, Imagerie Cerebrale et Handicaps Neurologiques UMR 825, Toulouse 2Universite de Toulouse, UPS, Imagerie Cerebrale et Handicaps Neurologiques UMR 825, CHU Purpan, Place du Dr Baylac, Toulouse Cedex 9 3LU 48 LERISM Laboratoire d’Etudes et de Recherche en Imagerie Spatiale et Medicale, UPS, Toulouse Cedex 4 4Universite de Toulouse, UPS, INPT, LAPLACE (Laboratoire Plasma et Conversion d'Energie), Toulouse Cedex 9 France


Introduction
Electroencephalography (EEG) and/or Event Related Potentials (ERP) are powerful noninvasive techniques which have broad clinical applications for epilepsy (Gloor et al., 1977;Hughes, 1989;Jaseja, 2009;Myatchin et al., 2009). It is also the case for psychiatric and developmental disorders (Pae et al., 2003;Ruchsow et al., 2003;Youn et al., 2003). There are developments in brain cognition research as for dyslexia (Horowitz-Kraus&Breznitz, 2008;Nuwer, 1998;Russeler et al., 2007), for visual treatment in face recognition (Chaby et al., 2003;George et al., 1996). In all these situations, specific brain areas are activated, and inverse techniques based on ERP treatment can help to estimate them. Techniques based on EEG/ERP are known to be incontestably inoffensive and cheap. This explains why they are often used and are still of great interest in medicine. The optimization of such medical tools, in research on brain cognition and/or as clinical tools, often requires knowledge of the intracerebral current sources. In EEG/ERP, this information can be obtained by solving of the socalled "inverse" problem consisting in finding the localization of the spatio-temporal intracerebral activity from scalp potential recordings. Various methods have been proposed in the EEG/ERP literature for computing this inverse problem.
Although scalp potentials were first recorded by Hans Berger in 1929(Berger, 1929, the first inverse problem approach was introduced by Cuffin et al. (Cuffin&Cohen, 1979) in both MEG and EEG, followed by Hämäläinen et al. in 1984(Hämäläinen&Ilmoniemi, 1984 in MEG. They later extended and increased the performance of the inverse approach applied to MEG (Hämäläinen&Ilmoniemi, 1994). The method was based on the Euclidean norm, which estimates the shortest vector solution in the source-current space (Hämäläinen&Ilmoniemi, 1994). This so-called Minimum Norm Estimate (MNE) is close to Tikhonov regularization (Tikhonov&Arsenin, 1977). However, the MNE solution is known to misreport actual deep sources as being in the outermost cortex (Pascual-Marqui, 1999;Pascual-Marqui et al., 2002). In order to compensate for the tendency of MNE to favour weak and surface sources, some authors have introduced a "weighting" matrix, calling this inverse method the Weighted Minimum Norm Estimate (WMNE) (Ding, 2009). Then, derived from this reasoning, many inverse methods have been used and/or improved specifically for EEG/ERP, modifying and/or reducing the solution space. Baillet introduced a priori to the solution which can be seen as a weighting matrix using a Bayesian probability, based on anatomical or functional knowledge (Baillet&Garnero, 1997). A Weighted Resolution Optimization (WROP), extending the Backus-Gilbert inverse method (Backus&Gilbert, 1968), has been developed (Grave de Peralta Menendez et al., 1997). The same technique has been modified, using biophysical and psychological a priori to the method called "Local Auto Regressive Average" (LAURA) (De Peralta-Menendez&Gonzales-Andino, 1998). Other authors have considered that restricting the potential solution to the cortical surface is sufficient to make the brain localization, and that the potential maps on the cortex surface must be significantly smooth, which has given rise to the inverse methods called "LOw Resolution brain Electromagnetic Tomography" (LORETA) (Pascual-Marqui et al., 1994), sLORETA (Pascual-Marqui, 2002 which is close to the "Variable Resolution Electrical Tomography" (VARETA) method (Bosch-Bayard et al., 2001). The above list of inverse methods is not exhaustive; a wide range of techniques exist for deriving inverse methods for use in EEG/ERP and new developments continue to be relevant today.
It should be noted that another type of inverse method has been developed at the same time. The main assumption is that the number of intra-cerebral current sources is limited (<10) and each source is punctual. Examples of such inverse methods are implemented in Brain Electric Source Analysis (BESA) (Scherg&Berg, 1991), using the so-called "simplex method" developed by Nelder and Mead (Nelder&Mead, 1965) and the Multiple Signal Classification (MUSIC) algorithm (Mosher et al., 1999). This type of method will not be discussed in our study, which only takes all cortex surfaces into account as possible locations of brain activity.
Inverse methods are numerous and cover many domains, especially in physics and medicine. Recent research can be found that uses the CGM for problems such as the determination of local boiling heat fluxes (Egger et al., 2009), the spatial distribution of Young's modulus (Fehrenbach et al., 2006), 3D elastic full-waveform seismic inversion (Epanomeritakis et al., 2008). Other applications can be found in other journals e.g., thermal diffusivity in plasma (Perez et al., 2008;Yang et al., 2008) and conductivity changes in impedance tomography (Zhao et al., 2007), proving, if it were necessary, the wide use of CGM in many different fields of application. Nevertheless, despite some attempts to use inverse methods such as the CGM in EEG/ERP, there is a lack of studies on the application of CGM to inverse problems in electroencephalography and/or event related potentials. Our contribution is to study the interest of applying CGM in EEG or ERP inverse problems.
In this article, the dependence of the reconstruction quality on the number of electrodes and the noise level are studied using CGM in numerical simulation. The main goal of this work is to evaluate the quality of intra-cerebral source reconstruction using CGM and to compare these results to the Cortical Imaging Technique (CIT). The model parameters and the CGM are described in Sec. 2. Then, in Sec. 3, the theoretical reconstruction of cortical potentials, as if they had been solved from experimentally recorded scalp potentials, are presented and www.intechopen.com discussed, considering various numbers of electrodes and noise levels. In Sec. 4, previous results are compared to those obtained by CIT, using the comparison tools MAG and RDM factors. The conclusions of this work are given in Sec. 5.

Head model
To localize brain activity from recorded scalp potentials in EEG/ERP, mathematical/ physical models that describe the geometrical and electrical properties of the head and the intra-cerebral current sources are needed. Generally, the head is described as a conductive volume with piecewise constant conductivity to represent the conductivity of each of its different parts. (Chauveau et al., 2004;He et al., 2002;Zhang et al., 2003). In our study (Figure 1), five compartments were used to construct the head model for the simulation, using the ICBM-152 (http://packages.bic.mni.mcgill.ca/tgz/) T1 template from Montreal Neurological Institute.
Resolution was 2 mm. Conductivities were those used in our previous study on CIT (Chauveau et al., 2008).  Fig. 1. The five tissues (white and grey matter, cerebrospinal fluid, skull and scalp), after segmentation of a realistic head geometry (e.g. The T1 Montreal head template).

Method
Determining scalp potentials from the simulation of intra-cerebral sources, called the forward problem, was an initial step towards the solution of the inverse problem, which aimed to find the sources at the origin of scalp potentials. Various numerical methods Darvas et al., 2006;Franceries et al., 2003) have been proposed in the literature for computing the forward problem, including finite difference (FDM) ( (Mattout, 2002;Vanrumste, 2001), boundary element (BEM) (Crouzeix, 2001;Kybic et al., 2005;Yvert et al., 1995) and finite element (FEM) (Darvas et al., 2006;Thevenet et al., 1991) methods, the last two being the most widely used. FEM with inclusion anisotropic conductivities have also been developed (Wolters et al., 2007). Although the simulations are usually time consuming, all give rise to numerical solutions and most of them are adequate to simulate brain activation. The Resistor Mesh Model (RMM) Franceries et al., 2003), close to Finite Volume Method (FVM) first proposed by Patankar (Patankar, 1980), gives very stable results and is easy to set up. The RMM is made of 2 mm size voxel elements. A sparse square symmetric admittance matrix Y describes the model. Each element represents a resistor, completely determined by its geometry and its conductivity (Franceries et al., 2003). The elements are assembled at the nodes of the model. The forward solution for a vector of currents I is a vector of potentials V so that I = Y x V. The resolution is obtained by using a numerical technique as Newton Raphson algorithm. It should be noted that, in some very special and simple cases (e.g. spherical models), an analytical solution is available but this is not the case for realistic head geometry (de Munck&Peters, 1993;Yvert et al., 1997;Zhou&van Oosterom, 1992)

Source configuration
In EEG/ERP, brain activation was first simulated by one or several dipolar current sources. The brain activation of each source was modelled by a current dipole, as introduced in 1953 by Plonsey (Plonsey&Barr, 1988). Other types of extended brain activity model have been proposed, e.g. ring extended sources to mimic the gamma frequency range EEG (Tallon-Baudry et al., 1999) .
In our study, we chose the current dipole model, which is widely used and well suited to the RMM. The intra-cerebral activation was simulated by four current dipoles, as used in a previous study on CIT (Chauveau et al., 2008), in order to make the comparison between CGM and CIT. A complex source configuration was used with 2 radial (RR = Radial Right, RL = Radial Left ) and 2 tangential ( TR = Tangential Right ant TL = Tangential Left) dipoles, placed on or close to the cortex surface. We chose these four dipoles because we wanted to test the inverse technique on two major points: radial or tangential dipoles (EEG being known to be most sensitive to radial), and symmetric dipoles (is the technique able to separate right and left activity?).

Forward solution
The RMM was applied to solve the forward problem with the previously described source configuration and a sample head model of five tissue compartments, using 2 mm voxels. The method computes potentials at all nodes (the RMN model contains 486,850 nodes and 1,413,720 elements) inside the head model and at the head surface where the electrodes are www.intechopen.com placed (i.e. on the patient's scalp surface). The use of a large matrix made solving this complete forward problem time consuming. In EEG/ERP, in order to reduce time of resolution and to minimize hard disk space, a lead field (LF) matrix, linking the electrode potentials and the currents at the cortex surface is constructed, using the Helmholtz reciprocity principle (Helmholtz, 1853). This new matrix is smaller, reducing calculation time, but the potentials are computed only at the electrode nodes of the scalp surface. The forward problem of computing scalp potential at electrode position (V e ) from a source configuration (I) thus becomes a reduced linear system as follows: .

Inverse problem
Generally, the inverse problem is solved by using the same matrix as the one for the numerical forward method (i.e. LF ), but using inversion. The LF matrix is not square and so cannot be inverted directly. Many methods exist to solve this ill-posed inverse problem, detailed mathematically by Tikhonov (Tikhonov&Arsenin, 1977). Depending on the physical problem, the matrix conditioning and the optimal inverse method have to be adapted. Up to now, CGM have not been applied to EEG/ERP and the most widely used inverse method is the pseudo-inverse matrix, e.g. the Moore-Penrose technique: where LF  and T LF are respectively the pseudo-inverse matrix of LF by Moore-Penrose and the transpose matrix of LF .
In real measurements, data are corrupted by noise and a regularization technique has to be used in the inversion procedure. Zero-order Tikhonov regularization permits this problem to be solved:  is a regularization factor depending on noise level, the optimal value of which is obtained at the angle of the associated L-curve (Carthy, 2003;Hansen, 2000;Tikhonov, 1963).
In EEG/ERP, the Cortical Imaging Technique (CIT) is one of the possible inverse methods, which limits the space of solutions for current dipoles to the cortex surface. This method has been described and evaluated (Chauveau et al., 2008;He et al., 2002) and it provided the comparison technique used in our study.

Conjugate gradient method (CGM)
CGM is an iterative technique. Other iterative techniques have been proposed (Gorodnitsky et al., 1995;Hansen, 1994;Ioannides et al., 1990). Ioannides proposes continuous probabilistic solutions to the biomagnetic inverse problem, very efficient for deep sources. Gorodnitsky describes a recursive weighted minimum norm algorithm (FOCUSS). Hansen has developed regularization tools for Matlab: he describes the iterative regularization methods, and presents CGM as a process which has some inherent regularization effect where the number of iterations plays the role of regularization parameter.
CGM was first designed for solving linear equations thanks to a square symmetric matrix. The application of CGM can be extended to rectangular non symmetric matrix as lead fields are, for the inverse solution. That is this way we use CGM in this study. The following equations do not need specific hypothesis on the properties of the linear matrix.
When using a gradient method (GM) in EEG/ERP, the inverse problem is replaced by an estimation problem in which the unknown source configuration k I is varied iteratively until the difference between the measured and calculated scalp potentials is as small as possible: .
The simple gradient method (Amari, 1977) is based on a local derivative function, in order to minimize the error. At each step of a gradient method, a trial set of values for the variable is used to generate a new set corresponding to a lower value of the error function.This was improved in the steepest gradient method (Curry, 1944), where the descent method takes the direction of the maximum gradient of the error function, which reduces the number of iterations. A further improvement is CGM, in which the previous (k) and the next (k+1) search directions are defined to be orthogonal in the residual associated error space, so that CGM explores a maximum of k R space. The CGM (Press et al., 1992) is an iterative method which computes: where k P is a vector of search direction at the k th iteration and  is a scalar of optimal step of descent obtained by finding the minimal argument of the objective function k F , defined by the norm of residual The solution k I is obtained at the k th iteration when the value of the chosen stopping criterion C of CGM is reached: C corresponds to a value, chosen by the user: it must be higher or equal to 0.01, from our experience and with our model. The root mean square of the relative error is compared to that value to stop the iterations.
In real conditions, data are corrupted by noise. In ERP/EEG protocols, noise vector No can be simply estimated on the pre-triggering time interval before events. Then the smaller criterion to reach becomes: . .
CGM does not need a priori conditions for solving the inverse problem in EEG/ERP, especially on the number of possible current sources in the cerebral volume. Moreover, CGM is faster than the classical Gradient Method because it needs less iteration to converge.
The main particularity of CGM is that the variation of the vector current obtained at each program loop is made orthogonal to the previous one. This permits to explore more quickly the space of solutions. In our application, CGM uses a lead field matrix, which is never inverted. The minimization is achieved on the square difference between measured and estimated electrode potentials. To stop the process, two methods are reported:  CGM performance is given for single and multiple dipoles in Table 1. It appears that the cortical potentials obtained by CGM underestimate the dipole amplitudes in comparison with potentials obtained by the forward solution. The worst result is observed for the dipole RR, which is correctly located on the cortex, but with a spread cortex area (Fig. 2) www.intechopen.com and with a local maximum much lower than for the forward solution, which explains the low value of MAG.

Effect of number of electrodes on CGM with the 4 dipoles
As EEG is only recorded at a limited number of electrodes, it is important to estimate the role of this number on the quality of the inverse solution. Figure 3 shows the cortical potentials obtained by CGM for 60 and 107 electrodes without noise. CGM was used to compute cortical potentials from the electrode potentials of the forward solution. As we can see on the figure, the CGM solution with 107 electrodes is more accurate than the solution obtained with 60 electrodes (tangential dipoles are better defined: red and blue areas are closer). We also observe, taking the potentials of the forward solution as a reference, that CGM with 60 and 107 electrodes underestimates the potentials at the cortical surface, especially for tangential dipoles. MAG reported in Table 1, and Table 2 and 3 (for noise 0%) confirms lower potential estimation at the cortex, whereas high RDM indicates mismatch on the shape or position.

CGM results with noise
A recent review on solving the inverse problem in EEG (Grech et al., 2008) presents the techniques in non-parametric and parametric methods, depending on the fixed number of dipoles (assumed a priori or not). No specific technique appears to give much better results than the others, and research in this field is continuing. For simulation studies, EEG noise must be taken into account, and Gaussian White Noise (GWN) is often used (Chauveau et al., 2008;He et al., 2002). We tested the CGM with 3 different GWN levels: 2%, 5% and 10% of the maximum electrode potential. Figure 4 shows the cortical potential distribution obtained by CGM with 60 and 107 electrodes for noise levels varying from 0% to 10% (criterion defined in equation 14).
Oscillations increase with the level of noise with 60 and 107 electrodes. So, the higher the noise level, the less correctly the cortical potential cartography is reconstructed. These results also show that oscillations in cortical potential distributions increase relatively faster when the noise level is higher than 5% with 107 electrodes, and 10% with 60 electrodes. Figure 5 shows the cortical potential distributions obtained with CGM, for different values of relative noise level and with 107 electrodes (criterion defined in equation 14), in comparison with the results of CIT. A Tikhonov regularization was used in CIT, but there are anyway some oscillations for high noise level. CGM presents oscillations when the criterion is too small compared to noise, but for each noise level it a correct estimation can be obtained.

CGM versus CIT
In real conditions, the noise level can be easily estimated on the pre-stimuli interval before the triggers used for ERP. Taking into account the noise level, the criterion is then limited to C noise (equation 15). Iterations are stopped when C noise is reached. Results are reported in figure 6.
Qualitative factors have been calculated by means of MAG and RDM (see appendix) at the electrodes, at the scalp surface and at the cortex surface for 60 electrodes (Table 2) and for 107 electrodes (Table 3).

Conclusion
This study by simulation has shown that CGM gives coherent results in the detection of simultaneous multiple dipoles (4 in our case). CGM solutions give satisfactory localization and estimation of cortical potentials even though the area of each dipole is increased. Symmetrical dipoles are well detected while tangential dipoles are more difficult to observe, as for any inverse technique in EEG/ERP.
We have shown that, without noise, CGM correctly localizes individual and simultaneous dipoles, with an underestimation of the cortical potentials. Moreover, the number of electrodes strongly conditions the quality of the solution obtained by CGM. So, without noise in the data, the higher the number of electrodes, the more accurate the dipole localization and the more correctly reconstructed the corresponding cortex potentials. So increasing the number of electrodes reduces the number of unknowns in the inverse problem in EEG/ERP. In consequence, cortical potentials are better evaluated.
With the addition of white Gaussian noise (WGN), this observation becomes partially true, because solutions obtained with high numbers of electrodes are less stable than those obtained with a smaller number when noise level increases. There is then an optimal number of electrodes for each simulated noise level. Solutions obtained by CIT and by CGM present oscillations which increase with the noise level. Cortical potential solutions of CGM are quite similar to the ones of CIT for a low noise level. When this level increases, CIT presents oscillation, still gives quite correct position for the sources but added potentials corrupt the result, while CGM presents solutions with less oscillation, but may be with less precision. The combination of CIT and CGM results permits to validate the source positions: CGM permit to clearly identify there are 4 sources in our case (2 radial dipoles and 2 tangential dipoles), and CIT permits to point out where they are.
It is then possible to use CGM as a complementary tool to solve inverse problems in EEG/ERP. The advantage of CGM is that there is no need for matrix inversion and there is not a prior in the number of current sources or in their propagation direction in the cerebral volume. This iterative method avoids having to invert huge rectangular matrices which are time and memory consuming when the spatial resolution of the model is ambitious. The estimation of noise permits to calculate a realistic stopping criterion to use, avoiding oscillations.

Appendix: Comparison tools MAG and RDM
MAG is an index for potential magnitude comparison between two series of equivalent data, and RDM estimates the variation of spatial distribution between the two series.
MAG and RDM are given by:

 
where Ci V is a series of computed potential data points obtained with a specific technique from electrode potentials ( CIT or CGM), Fi V is the forward solution for the same points and n is the number of chosen points.