 Review
 Open Access
Review on solving the forward problem in EEG source analysis
 Hans Hallez^{1}Email author,
 Bart Vanrumste^{2, 3}Email author,
 Roberta Grech^{4},
 Joseph Muscat^{6},
 Wim De Clercq^{2},
 Anneleen Vergult^{2},
 Yves D'Asseler^{1},
 Kenneth P Camilleri^{5},
 Simon G Fabri^{5},
 Sabine Van Huffel^{2} and
 Ignace Lemahieu^{1}
https://doi.org/10.1186/17430003446
© Hallez et al; licensee BioMed Central Ltd. 2007
Received: 05 January 2007
Accepted: 30 November 2007
Published: 30 November 2007
Abstract
Background
The aim of electroencephalogram (EEG) source localization is to find the brain areas responsible for EEG waves of interest. It consists of solving forward and inverse problems. The forward problem is solved by starting from a given electrical source and calculating the potentials at the electrodes. These evaluations are necessary to solve the inverse problem which is defined as finding brain sources which are responsible for the measured potentials at the EEG electrodes.
Methods
While other reviews give an extensive summary of the both forward and inverse problem, this review article focuses on different aspects of solving the forward problem and it is intended for newcomers in this research field.
Results
It starts with focusing on the generators of the EEG: the postsynaptic potentials in the apical dendrites of pyramidal neurons. These cells generate an extracellular current which can be modeled by Poisson's differential equation, and Neumann and Dirichlet boundary conditions. The compartments in which these currents flow can be anisotropic (e.g. skull and white matter). In a threeshell spherical head model an analytical expression exists to solve the forward problem. During the last two decades researchers have tried to solve Poisson's equation in a realistically shaped head model obtained from 3D medical images, which requires numerical methods. The following methods are compared with each other: the boundary element method (BEM), the finite element method (FEM) and the finite difference method (FDM). In the last two methods anisotropic conducting compartments can conveniently be introduced. Then the focus will be set on the use of reciprocity in EEG source localization. It is introduced to speed up the forward calculations which are here performed for each electrode position rather than for each dipole position. Solving Poisson's equation utilizing FEM and FDM corresponds to solving a large sparse linear system. Iterative methods are required to solve these sparse linear systems. The following iterative methods are discussed: successive overrelaxation, conjugate gradients method and algebraic multigrid method.
Conclusion
Solving the forward problem has been well documented in the past decades. In the past simplified spherical head models are used, whereas nowadays a combination of imaging modalities are used to accurately describe the geometry of the head model. Efforts have been done on realistically describing the shape of the head model, as well as the heterogenity of the tissue types and realistically determining the conductivity. However, the determination and validation of the in vivo conductivity values is still an important topic in this field. In addition, more studies have to be done on the influence of all the parameters of the head model and of the numerical techniques on the solution of the forward problem.
Keywords
Introduction
Since the 1930s electrical activity of the brain has been measured by surface electrodes connected to the scalp [1]. Potential differences between these electrodes were then plotted as a function of time in a socalled electroencephalogram (EEG). The information extracted from these brain waves was, and still is instrumental in the diagnoses of neurological diseases [2], mainly epilepsy. Since the 1960s the EEG was also used to measure eventrelated potentials (ERPs). Here brain waves were triggered by a stimulus. These stimuli could be of visual, auditory and somatosensory nature. Different ERP protocols are now routinely used in a clinical neurophysiology lab. Researchers nowadays are still searching for new ERP protocols which may be able to distinguish between ERPs of patients with a certain condition and ERPs of normal subjects. This could be instrumental in disorders, such as psychiatric and developmental disorders, where there is often a lack of biological objective measures.
During the last two decades, increasing computational power has given researchers the tools to go a step further and try to find the underlying sources which generate the EEG. This activity is called EEG source localization. It consists of solving a forward and inverse problem. Solving the forward problem starts from a given electrical source configuration representing active neurons in the head. Then the potentials at the electrodes are calculated for this configuration. The inverse problem attempts to find the electrical source which generates a measured EEG. By solving the inverse problem, repeated solutions of the forward problem for different source configurations are needed. A review on solving the inverse problem is given in [3].
In this review article several aspects of solving the forward problem in EEG source localization will be discussed. It is intended for researchers new in the field to get insight in the stateoftheart techniques to solve the forward problem in EEG source analysis. It also provides an extensive list of references to the work of other researchers.
First, the physical context of EEG source localization will be elaborated on and then the derivation of Poisson's equation with its boundary conditions. An analytical expression is then given for a threeshell spherical head model. Along with realistic head models, obtained from medical images, numerical methods are then introduced that are necessary to solve the forward problem. Several numerical techniques, the Boundary Element Method (BEM), the Finite Element Method (FEM) and the Finite Difference Method (FDM), will be discussed. Also anisotropic conductivities which can be found in the white matter compartment and skull, will be handled.
The reciprocity theorem used to speed up the calculations, is discussed. The electric field that results at the dipole location within the brain due to current injection and withdrawal at the surface electrode sites is first calculated. The forward transfercoefficients are obtained from the scalar product of this electric field and the dipole moment. Calculations are thus performed for each electrode position rather than for each dipole position. This speeds up the time necessary to do the forward calculations since the number of electrodes is much smaller than the number of dipoles that need to be calculated.
The number of unknowns in the FEM and FDM can easily exceed the million and thus lead to large but sparse linear systems. As the number of unknowns is too large to solve the system in a direct manner, iterative solvers need to be used. Some popular iterative solvers are discussed such as successive overrelaxation (SOR), conjugate gradient method (CGM) and algebraic multigrid methods (AMG).
The physics of EEG
In this section the physiology of the EEG will be shortly described. In our opinion, it is important to know the underlying mechanisms of the EEG. Moreover, forward modeling also involves a good model for the generators of the EEG. The mechanisms of the neuronal actionpotentials, excitatory postsynaptic potentials and inhibitory postsynaptic potentials are very complex. In this section we want to give a very comprehensive overview of the underlying neurophysiology.
Neurophysiology
The brain consists of about 10^{10} nerve cells or neurons. The shape and size of the neurons vary but they all possess the same anatomical subdivision. The soma or cell body contains the nucleus of the cell. The dendrites, arising from the soma and repeatedly branching, are specialized in receiving inputs from other nerve cells. Via the axon, impulses are sent to other neurons. The axon's end is divided into branches which form synapses with other neurons. The synapse is a specialized interface between two nerve cells. The synapse consists of a cleft between a presynaptic and postsynaptic neuron. At the end of the branches originating from the axon, the presynaptic neuron contains small rounded swellings which contain the neurotransmitter substance. Further readings on the anatomy of the brain can be found in [4] and [5].
One neuron generates a small amount of electrical activity. This small amount cannot be picked up by surface electrodes, as it is overwhelmed by other electrical activity from neighbouring neuron groups. When a large group of neurons is simultaneously active, the electrical activity is large enough to be picked up by the electrodes at the surface and thus generating the EEG. The electrical activity can be modeled as a current dipole. The current flow causes an electric field and also a potential field inside the human head. The electric field and potential field spreads to the surface of the head and an electrode at a certain point can measure the potential [2].
At rest the intracellular environment of a neuron is negatively polarized at approximately 70 mV compared with the extracellular environment. The potential difference is due to an unequal distribution of Na^{+}, K^{+} and Cl^{} ions across the cell membrane. This unequal distribution is maintained by the Na^{+} and K^{+} ion pumps located in the cell membrane. The GoldmanHodgkinKatz equation describes this resting potential and this potential has been verified by experimental results [2, 6, 7].
The neuron's task is to process and transmit signals. This is done by an alternating chain of electrical and chemical signals. Active neurons secrete a neurotransmitter, which is a chemical substance, at the synaptical side. The synapses are mainly localized at the dendrites and the cell body of the postsynaptic cell. A postsynaptic neuron has a large number of receptors on its membrane that are sensitive for this neurotransmitter. The neurotransmitter in contact with the receptors changes the permeability of the membrane for charged ions. Two kinds of neurotransmitters exist. On the one hand there is a neurotransmitter which lets signals proliferate. These molecules cause an influx of positive ions. Hence depolarization of the intracellular space takes place. A depolarization means that the potential difference between the intra and extracellular environment decreases. Instead of 70 mV the potential difference becomes 40 mV. This depolarization is also called an excitatory postsynaptic potential (EPSP). On the other hand there are neurotransmitters that stop the proliferation of signals. These molecules will cause an outflow of positive ions. Hence a hyperpolarization can be detected in the intracellular volume. A hyperpolarization means that the potential difference between the intra and extracellular environment increases. This potential change is also called an inhibitory postsynaptic potential (IPSP). There are a large number of synapses from different presynaptic neurons in contact with one postsynaptic neuron. At the cell body all the EPSP and IPSP signals are integrated. When a net depolarization of the intracellular compartment at the cell body reaches a certain threshold, an action potential is generated. An action potential then propagates along the axon to other neurons [2, 6, 7].
The generators of the EEG
The electrodes used in scalp EEG are large and remote. They only detect summed activities of a large number of neurons which are synchronously electrically active. The action potentials can be large in amplitude (70–110 mV) but they have a small time course (0.3 ms). A synchronous firing of action potentials of neighboring neurons is unlikely. The postsynaptic potentials are the generators of the extracellular potential field which can be recorded with an EEG. Their time course is larger (10–20 ms). This enables summed activity of neighboring neurons. However their amplitude is smaller (0.1–10 mV) [3, 8].
Apart from having more or less synchronous activity, the neurons need to be regularly arranged to have a measurable scalp EEG signal. The spatial properties of the neurons must be so that they amplify each other's extracellular potential fields. The neighboring pyramidal cells are organized so that the axes of their dendrite tree are parallel with each other and normal to the cortical surface. Hence, these cells are suggested to be the generators of the EEG.
A migration of positively charged ions from the cell body and the basal dendrites to the apical dendrite occurs, which is illustrated in figure 2(a) with current lines. This configuration generates extracellular potentials. Other membrane activities start to compensate for the massive intrusion of the positively charged ions at the apical dendrite, however these mechanisms are beyond the scope of this work and can be found elsewhere [2, 9, 10].
A simplified equivalent electric circuit is presented in figure 2(b) to illustrate the initial activity of an EPSP. At rest, the potential difference between the intra and extracellular compartments can be represented by charged capacitors. One capacitor models the potential difference at the apical dendrites side while a second capacitor models the potential difference at the cell body and basal dendrite side. The potential difference over the capacitors is 60 mV. The neurotransmitter causes a massive intrusion of positively charged ions at the postsynaptic membrane at the apical dendrite side. In the equivalent circuit, this is modeled by a switch that is closed. The capacitor at the cell body side discharges causing a current flow over the extracellular resistor R_{ e }and the intracellular resistor R_{ i }. The repolarization of the cell membrane at the apical side or the initiation of the action potential are not modeled with this simple equivalent electrical circuit.
The capacitors and the switch, in figure 2(b), represent a model of the electrical source at the initial phase of the depolarization of the neuron. They could also be replaced by a time dependent current source, however this representation is not ideal. The capacitor representation, for the initial phase of depolarization, fits closer the occurring physical phenomena. The impedance of the tissue in the human head has, for the frequencies contained in the EEG, no capacitive nor inductive component and is hence pure resistive. More advanced equivalent electrical circuits can be found elsewhere [10]. The fact that a current flows through the extracellular resistor indicates that potential differences in the extracellular space can be measured.
A simplified electrical model for this active cell consists of two current monopoles: a current sink at the apical dendrite side which removes positively charged ions from the extracellular environment, and a current source at the cell body side which injects positively charged ions in the extracellular environment. The extracellular resistance R_{ e }can be decomposed in the volume conductor model in which the active neuron is embedded, as illustrated in figure 2(c). For further reading on the generation of the EEG one can refer to [11] and [9].
Poisson's equation, boundary conditions and dipoles
In the previous sections we saw that the generators of the EEG are the synaptic potentials along the apical dendrites of the pyramidal cells of the grey matter cortex. It is important to notice that the EEG reflects the electrical activity of a subgroup of neurons, especially pyramidal neuron cells, where the apical dendrite is systematically oriented orthogonal to the brain surface. Certain types of neurons are not systematically oriented orthogonal to the brain surface. Therefore, the potential fields of the synaptic currents at different dendrites of neurons van cancel each other out. In that case the neuronal activity is not visible at the surface. Moreover, that actionpotentials, propagating along the axons, have no influence on the EEG. Their short timespan (2 ms) make the chance of generating simultaneous actionpotentials very small [6, 12]. In this section, a mathematical approach on the generation of the forward problem is given.
Quasistatic conditions
It is shown in [13] that no charge can be piled up in the conducting extracellular volume for the frequency range of the signals measured in the EEG. At one moment in time all the fields are triggered by the active electric source. Hence, no time delay effects are introduced. All fields and currents behave as if they were stationary at each instance. These conditions are also called quasistatic conditions. They are not static because the neural activity changes with time. But the changes are slow compared to the propagation effects.
Applying the divergence operator to the current density
Poisson's equation gives a relationship between the potentials at any position in a volume conductor and the applied current sources. The mathematical derivation of Poisson's equation via Maxwell's equations, can be found in various textbooks on electromagnetism [6, 10, 14]. Poisson's equation is derived with the divergence operator. In this way the emphasis is, in our opinion, more on the physical aspect of the problem. Furthermore, the concepts introduced in [10, 14], such as current source and current sink, are used when applying the divergence operator.
Definition
The integral over a closed surface ∂G represents a flux or a current. This integral is positive when a net current leaves the volume G and is negative when a net current enters the volume G. The vector d S for a surface element of ∂G with area dS and outward normal e_{ n }, can also be written as e_{ n }dS. The unit of ∇·J is A/m^{3} and is often called the current source density which in [15] is symbolized with I_{ m }. Generally one can write:
∇·J = I_{ m }.
Applying the divergence operator to the extracellular current density
First a small volume in the extracellular space, which encloses a current source and current sink, is investigated. The current flowing into the infinitely small volume, must be equal to the current leaving that volume. This is due to the fact that no charge can be piled up in the extracellular space. The surface integral of equation (1) is then zero, hence ∇·J = 0.
In the second case a volume enclosed by the current sink with position parameters r_{1}(x_{1}, y_{1}, z_{1}) is assumed. The current sink represents the removal of positively charged ions at the apical dendrite of the pyramidal cell. The integral of equation (1) remains equal to I while the volume G in the denominator becomes infinitesimally small. This gives a singularity for the current source density. This singularity can be written as a delta function: Iδ(r  r_{1}). The negative sign indicates that current is removed from the extracellular volume. The delta function indicates that current is removed at one point in space.
Uniting the three cases given above, one obtains:
∇·J = Iδ(r  r_{2})  Iδ(r  r_{1}).
Ohm's law, the potential field and anisotropic/isotropic conductivities
The relationship between the current density J in A/m^{2} and the electric field E in V/m is given by Ohm's law:
J = σ E
and with units A/(Vm) = S/m. There are tissues in the human head that have an anisotropic conductivity. This means that the conductivity is not equal in every direction and that the electric field can induce a current density component perpendicular to it with the appropriate σ in equation (4).
White matter consists of different nerve bundles (groups of axons) connecting cortical grey matter (mainly dendrites and cell bodies). The nerve bundles consist of nerve fibres or axons (see figure 4(b)). Water and ionized particles can move more easily along the nerve bundle than perpendicular to the nerve bundle. Therefore, the conductivity along the nerve bundle is measured to be 9 times higher than perpendicular to it [19, 20]. The nerve bundle direction can be estimated by a recent magnetic resonance technique: diffusion tensor magnetic resonance imaging (DTMRI) [21]. This technique provides directional information on the diffusion of water. It is assumed that the conductivity is the highest in the direction in which the water diffuses most easily [22]. Authors [23–25] have showed that anisotropic conducting compartments should be incorporated in volume conductor models of the head whenever possible.
The reference values of the absolute and relative conductivity of the compartments incorporated in the human head.
compartments  Geddes & Baker (1967)  Oostendorp (2000)  Gonçalves (2003)  Guttierrez (2004)  Lai (2005) 

scalp  0.43  0.22  0.33  0.749  0.33 
skull  0.006 – 0.015  0.015  0.0081  0.012  0.0132 
cerebrospinal fluid        1.79   
brain  0.12 – 0.48  0.22  0.33  0.313  0.33 
σ_{ scalp }/σ_{ skull }  80  15  20–50  26  25 
The skull conductivity has been subject to debate among researchers. In vivo measurements are very different from in vitro measurements. On top of that, the measurements are very patient specific. In [27], it was stated that the skull conductivity has a large influence on the forward problem.
It was believed that the conductivity ratio between skull and soft tissue (scalp and brain) was on average 80 [20]. Oostendorp et al. used a technique with realistic head models by which they passed a small current by means of 2 electrodes placed on the scalp. A potential distribution is then generated on the scalp. Because the potential values and the current source and sink are known, only the conductivities are unknown in the head model and equation (4) can be solved toward σ. Using this technique they could estimate the skulltosoft tissue conductivity ratio to be 15 instead of 80 [28]. At the same time, Ferree et al. did a similar study using spherical head models. Here, skulltosoft tissue conductivity was calculated as 25. It was shown in [29] that using a ratio of 80 instead of 16, could yield EEG source localization errors of an average of 3 cm up to 5 cm.
One can repeat the previous experiment for a lot of different electrode pairs and an image of the conductivity can be obtained. This technique is called electromagnetic impedance tomography or EIT. In short, EIT is an inverse problem, by which the conductivities are estimated. Using this technique, the skulltosoft tissue conductivity ratio was estimated to be around 20–25 [30, 31]. However in [30], it was shown that the skulltosoft tissue ratio could differ from patient to patient with a factor 2.4. In [32], maximum likelihood and maximum a posteriori techniques are used to simultaneously estimate the different conductivities. There they estimated the skulltosoft tissue ratio to be 26.
Another study came to similar results using a different technique. In Lai et al., the authors used intracranial and scalp electrodes to get an estimation of the skulltosoft tissue ratio conductivity. From the scalp measures they estimated the cortical activity by means of a cortical imaging technique. The conductivity ratio was adjusted so that the intracranial measurements were consistent with the result of the imaging from the scalp technique. They resulted in a ratio of 25 with a standard deviation of 7. One has to note however that the study was performed on pediatric patients which had the age of between 8 and 12. Their skull tissue normally contains a larger amount of ions and water and so may have a higher conductivity than the adults calcified cranial bones [33]. In a more experimental setting, the authors of [34] performed conductivity measures on the skull itself in patients undergoing epilepsy surgery. Here the authors estimated the skull conductivity to be between 0.032 and 0.080 S/m, which comes down to a softtissue to skull conductivity of 10 to 40.
Poisson's equation
The scalar potential field V, having volt as unit, is now introduced. This is possible due to Faraday's law being zero under quasistatic conditions (∇ × E = 0) [35]. The link between the potential field and the electric field is given utilizing the gradient operator,
E = ∇V.
The vector ∇V at a point gives the direction in which the scalar field V, having volt as its unit, most rapidly increases. The minus sign in equation (6) indicates that the electric field is oriented from an area with a high potential to an area with a low potential. Figure 3 also illustrates some equipotential lines generated by a current source and current a sink.
When equation (2), equation (4) and equation (6) are combined, Poisson's differential equation is obtained in general form:
∇·(σ∇(V)) = I_{ m }.
For the problem at hand, equation (3), equation (4) and equation (6) are combined yielding:
∇·(σ∇(V)) = Iδ(r  r_{2}) + Iδ(r  r_{1}).
The potentials V are calculated with equations (8), (9) or (10) for a given current source density I_{ m }, in a volume conductor model, e.g. in our application, the human head. Compartments in which all conductivities are equal, are called homogeneous conducting compartments.
Boundary conditions
where e_{ n }is the normal component on the interface.
Equations (11) and (12) are called the Neumann boundary condition and the homogeneous Neumann boundary condition, respectively.
The second boundary condition only holds for interfaces not connected with air. By crossing the interface the potential cannot have discontinuities,
V_{1} = V_{2}.
This equation represents the Dirichlet boundary condition.
The current dipole
The dipole moment d is defined by a unit vector e_{ d }(which is directed from the current sink to the current source) and a magnitude given by d = d = I·p, with p the distance between the two monopoles. Hence one can write:
d = I·p e_{ d }.
It is often so that a dipole is decomposed in three dipoles located at the same position of the original dipole and each oriented along one of the Cartesian axes. The magnitude of each of these dipoles is equal to the orthogonal projection on the respective axis as illustrated in figure 6(b). one can write:
d = d_{ x }e_{ x }+ d_{ y }e_{ y }+ d_{ z }e_{ z },
with e_{ x }, e_{ y }and e_{ z }being the unit vectors along the three axes. Furthermore, d_{ x }, d_{ y }and d_{ z }are often called the dipole components. Notice that Poisson's equation (8) is linear. Due to a dipole at a position r_{ dip }and dipole moment d, a potential V at an arbitrary scalp measurement point r can be decomposed in:
V(r, r_{ dip }, d) = d_{ x }V(r, r_{ dip }, e_{ x }) + d_{ y }V(r, r_{ dip }, e_{ y }) + d_{ z }V(r, r_{ dip }, e_{ z }).
A large group of pyramidal cells need to be more or less synchronously active in a cortical patch to have a measurable EEG signal. All these cells are furthermore oriented with their longitudinal axis orthogonal to the cortical surface. Due to this arrangement the superposition of the individual electrical activity of the neurons results in an amplification of the potential distribution. A large group of electrically active pyramidal cells in a small patch of cortex can be represented as one equivalent dipole on macroscopic level [36, 37]. It is very difficult to estimate the extent of the active area of the cortex as the potential distribution on the scalp is almost identical to that of an equivalent dipole [38].
General algebraic formulation of the forward problem
In symbolic terms, the EEG forward problem is that of finding, in a reasonable time, the scalp potential g(r, r_{ dip }, d) at an electrode positioned on the scalp at r due to a single dipole with dipole moment d = d e_{ d }(with magnitude d and orientation e_{ d }), positioned at r_{ dip }. This amounts to solving Poisson's equation to find the potentials V(r) on the scalp for different configurations of r_{ dip }and d. For multiple dipole sources, the electrode potential would be $V(r)={\displaystyle \sum _{i}g(r,{r}_{di{p}_{i}},{d}_{i})=}{\displaystyle \sum _{i}g(r,{r}_{di{p}_{i}},{e}_{{d}_{i}}){d}_{i}}$. In practice, one calculates a potential between an electrode and a reference (which can be another electrode or an average reference).
where i = 1,...,p and j = 1,...,N. Here V is a column vector.
where V is now the matrix of data measurements, G is the gain matrix and D is the matrix of dipole magnitudes at different time instants.
More generally, a noise or perturbation matrix n is added,
V = GD + n.
In general for simulations and to measure noise sensitivity, noise distribution is a gaussian distribution with zero mean and variable standard deviation. However in reality, the noise is coloured and the distribution of the frequency depends on a lot of factors: patient, measurement setup, pathology,....
A general multipole expansion of the source model
Solving the inverse problem using multiple dipole model requires the estimation of a large number of parameters, 6 for each dipole. Given the use of a limited number of EEG electrodes, the problem becomes underdetermined. In this case, regularization techniques have to be applied, but this leads to oversmoothed source estimations. On the other hand, the use of a limited number of dipoles (one, two or three) leads to very simplified sources, which are very often ambiguous and cause errors due to simplified modelling. The dipole model as a source is a good model for focal brain activity.
A multipole expansion is an alternative (first introduced by [39]), which is based on a spherical harmonic expansion of the volume source, which is not necessarily focal. It provides the added model flexibility needed to account for a wide range of physiologically plausible sources, while at the same time keeping the number of estimation parameters sufficiently low. In fact, The zerothorder and firstorder terms in the expansion are called the monopole and dipole moment, respectively. A quadrupole is a higher order term and is generated by two equal and oppositely oriented dipoles whose moments tend to infinity as they are brought infinitesimally close to each other. An octapole consists of two quadrupoles brought infinitesimally close to each other and so on. It can be shown that if the volume G containing the active sources I_{ sv }(r ') is limited in extent, the solution to Poisson's equation for the potential V may be expanded in terms of a multipole series:
V = V_{ monopole }+ V_{ dipole }+ V_{ quadrupole }+ V_{ octapole }+ V_{ hexadecpole }+ ...
where V_{ quadrupole }is the potential field caused by the quadrupole. In practice, a truncated multipole series is used up to a quadrupole, because the contribution to the electrode potentials by a octapole or higher order sources rapidly decreases when the distance between electrode and source is increasing. The use of quadrupoles can sound plausible in the following case: A traveling action potential causes a depolarization wave through the axon, followed by a repolarization wave. These two phenomenon produce two opposite oriented dipoles very close to each other [40]. In sulci, pyramidal cells are oriented toward each other, which makes the use of quadrupole also reasonable. However, the skull causes a strong attenuation of the electrical field created by the source. Therefore, even a quadrupole has low contribution to the electrode potentials of the EEG, created by the volume current in the extracellular region. In EEG and ECG multiple dipoles of dipole layers are preferred over a multipole. Multipoles are popular in magnetoencephalographic (MEG) source localization, because of its low sensitivity to the skull conductivity [6, 10, 41, 42].
Solving the forward problem
Dipole field in an infinite homogeneous isotropic conductor
Equation (20) shows that a dipole field attenuates with 1/r^{2}. It is significant to remark that V, from equation (19), added with an arbitrary constant, is also a solution of Poisson's equation. A reference potential must be chosen. One can choose to set one electrode to zero or one can opt for average referenced potentials. The latter result in electrode potentials that have a zero mean.
The spherical head model
Where:
d_{ r }is the radial component (3 × 1vector in meters),
d_{ t }is the tangential component (3 × 1vector in meters),
R is the radius of the outer shell (meters),
S is the conductivity of scalp and brain tissue (Siemens/meter),
X is the ratio between the skull and soft tissue conductivity (unitless),
b is the relative distance of the dipole from the centre (unitless),
θ is the polar angle of the surface point see figure 8 (radians),
P_{ i }(·) is the Legendre polynomial,
${P}_{i}^{1}(\cdot )$ is the associated Legendre polynomial,
i is an index,
i_{1} equals 2i + 1,
r_{1} is the radius of the inner shell (in meters),
r_{2} is the radius of the middle shell (in meters),
f_{1} equals r_{1}/R (unitless) and
f_{2} equals r_{2}/R (unitless).
Equation (21) gives the scalp potentials generated by a dipole located on the zaxis, with zero dipole moment in the y direction. To find the scalp potentials generated by an arbitrary dipole, the coordinate system has to be rotated accordingly. Typical radii of the outer boundaries of the brain, skull and scalp compartments are equal to 8 cm, 8.5 cm and 9.2 cm, respectively [46]. An illustration of a typical spherical head model is shown in figure 8. These radii can be altered to fit a sphere more to the human head. The infinite series of equation (21) is often truncated. If the first 40 terms are used, the maximum scalp potential obtained with the truncated series, deviates less than 0.1% from the case where 100 terms are applied, for dipoles with a radial position smaller than 95% of the maximum brain radius.
There are also semianalytical solutions available for layered spheroidal anisotropic volume conductors [47–49]. Here the conductivity in the tangential direction can be chosen differently than in the radial direction of the sphere. Analytic solutions also exist for prolate and oblate spheroids or eccentric spheres [50–52].
Variants of the threeshell spherical head model, such as the Berg approximation [53], in which a singlesphere model is used to approximate a three (or four) layer sphere, have also been used to improve further the computational efficiency of multilayer spherical models.
Recently however, it is becoming more apparent that the actual geometry of the head [54–56] together with the varying thickness and curvatures of the skull [57, 58], affects the solutions appreciably. Socalled real head models are becoming much more common in the literature, in conjunction with either boundaryelement, finiteelement, or finitedifference methods. However, the computational requirements for a realistic head model are higher than that for a multilayer sphere.
An approach which is situated between the spherical head model approaches and realistic ones is the sensorfitted sphere approach [59]. Here a multilayer sphere is fitted to each sensor located on the surface of a realistic head model.
The boundary element method
The boundary element method (BEM) is a numerical technique for calculating the surface potentials generated by current sources located in a piecewise homogeneous volume conductor. Although it restricts us to use only isotropic conductivities, it is still widely used because of its low computational needs. The method originated in the field of electrocardiography in the late sixties and made its entrance in the field of EEG source localization in the late eighties [60]. As the name implies, this method is capable of providing a solution to a volume problem by calculating the potential values at the interfaces and boundary of the volume induced by a given current source (e.g. a dipole). The interfaces separate regions of differing conductivity within the volume, while the boundary is the outer surface seperating the nonconducting air with the conducting volume.
In practice, a head model is built from surfaces, each encapsulating a particular tissue. Typically, head models consist of 3 surfaces: brainskull interface, skullscalp interface and the outer surface. The regions between the interfaces are assumed to be homogeneous and isotropic conducting. To obtain a solution in such a piecewise homogenous volume, each interface is tesselated with small boundary elements.
where σ_{0} corresponds to the medium in which the dipole source is located (the brain compartment) and V_{0}(r) is the potential at r for an infinite medium with conductivity σ_{0} as in equation (19). ${\sigma}_{j}^{}$ and ${\sigma}_{j}^{+}$ are the conductivities of the, respectively, inner and outer compartments divided by the interface S_{ j }. d S is a vector oriented orthogonal to a surface element and d S the area of that surface element.
This equation can be transformed into a set of linear equations:
V = BV + V_{0},
where V and V_{0} are column vectors denoting at every node the wanted potential value and the potential value in an infinite homogeneous medium due to a source, respectively. B is a matrix generated from the integrals, which depends on the geometry of the surfaces and the conductivities of each region.
where r_{ i }, r_{ j }, r_{ k }are the nodes of the triangle and the triple scalar product is defined as [r_{ i }r_{ j }r_{ k }] = det(r_{ i }, r_{ j }, r_{ k }). The notation Δ_{i(jk) }is used to indicate any triangle for which one vertex is defined by the vector r_{ i }, the remaining two vertices denoted as r_{ j }and r_{ k }. The function h_{ i }(r) attains a value of unity at the ith vertex and drops linearly to zeros at the opposite edge of all triangles to which r_{ i }is a vertex. In this case, the collocation points are the vertices of the elements [66]. The approaches can be expanded into higherorder elements [67]. Gençer and Tanzer investigated quadratic and cubic element types and concluded that these gave superior results to models with linear elements [68].
where e is a vector with all its N (the total number of unknowns) components equal to one. The deflated equation
V = CV + V_{0},
possesses a unique solution which is also a solution to the orignal equation (27). If I denotes the N × N identity matrix and A represents I  C then
V = A^{1}V_{0}.
This equation can be solved using direct or iterative solvers. Direct solvers are especially usefull when the matrix A is relatively small because of a coarse grid. If one wants to use a fine grid, then iterative methods should be used. The use of multiple deflations during the iterations can significantly increase the convergence time to the solution of equation (31) [69].
A typical head model for solving the forward problem involves 3 layers: the brain, the skull and the scalp. The conductivity of the skull is lower than the conductivity of brain and scalp. If β is defined as the ratio of the skull conductivity to the brain conductivity Meijs et al. showed that an accurate solution of equation (23) is difficult to obtain for small β (β < 0.1). The large difference between the conductivities will cause an amplification of the numerical errors in the calculation. To solve this problem, the Isolated Problem Approach (IPA) can be used (also called Isolated Skull Approach), which was introduced by Hämäläinen and Sarvas [70]. Assume the labeling of the compartments as C_{ scalp }, C_{ skull }and C_{ brain }and S_{ scalp }as the outer interface, S_{ skull }as the interface between C_{ skull }and C_{ brain }and S_{ brain }as the interface between C_{ brain }and C_{ skull }. The IPA uses the following decomposition of the potential values:
V(r) = V'(r) + V''(r)
where ${V}_{0}^{1}$, ${V}_{0}^{2}$ and ${V}_{0}^{3}$ are the potentials at respectively the brainskull surface, the skullscalp surface and the outer surface. This imposes that ${{{V}^{\prime}}^{\prime}}_{0}^{3}$ has to be calculated. This can be done by solving the potentials at S_{ brain }with the scalp and skull compartments omitted. The increase in accuracy comes at a small cost of computational speed. A weighted IPA approach was developed by Fuchs et al. [71]. The IPA approach was extended to multisphere models by Gençer and AkahnAcar [72]. The calculation of the forward problem involves every node on the mesh, making it very computation intesive. Accelerated BEM computes the node potentials on a small subset of nodes corresponding to the electrode positions [73].
To improve the localization accuracy, one can locally refine the mesh. Yvert et al. showed that if the dipole is at 2 cm below the surface, a mesh of 0.5 triangles/cm^{2} is needed to have acceptable results. However, for shallow dipoles (between 2 mm and 20 mm below the brain surface) a mesh density of 2–6 triangles per cm^{2} is needed to obtain comparable results. Of course, the area in the mesh that has to be refined, has to be defined.
A main disadvantage using BEM in the EEG forward problem is that in all aforementioned implementations the precision drops when the distance of the source to one of the surfaces becomes comparable to the size of the triangles in the mesh. Kybic et al. presented a new framework based on a theorem that characterizes harmonic functions defined on the complement of a bounded smooth surface [74]. Using this framework, they proposed a symmetric formulation. The main benefit of this approach is that the error increases much less dramatically when the current sources approach a surface where the conductivity is discontinuous. In another paper by the same authors, a fast multipole acceleration was used to overcome the complexity of the symmetric formulation [75]. A recent article of the same authors demonstrates that the framework allows the use of more realistical head models, which don't have to be nested. In nested head models, an inner interface is completely enveloped by an outer interface. Nonnested compartments are compartments that are not part of the brain, but part of the head (such as eyes, sinuses,...) [76].
The finite element method
If (v, w) = ∫_{ G }v(x, y, z)w(x, y, z)dG and a(u, v) = (∇v, σ∇u), this can be written as:
a(V, φ) = (I_{ m }, φ)
where ${\left\{{\phi}_{i}\right\}}_{i=1}^{n}$ denotes a set of test functions also called basis functions. They have a local support, i.e. the area in which they are nonzero is limited to adjacent elements. Moreover, the basis functions span a space of piecewise polynomial functions.
Due to the local support of the basis function, each equation consists only of a linear combination of V_{ i }'s and its adjacent computational points. Hence the system A ∈ ℝ^{n×n}, A_{ ij }= a(φ_{ i }, φ_{ j }) is sparse. In matrix notation one can obtain:
A·V = I,
with I ∈ ℝ^{n×1 }being the column vector of the source terms obtained by the right hand side of equation (41).
An important consideration in finiteelement methods is how to represent a dipole source in the model.
• The obvious direct method is to represent a dipole using a pair of fixed voltage conditions of opposite polarity applied to two adjacent nodes [78].
• Another method is to embed a dipole source in the element basis functions. When the dipole lies along the edge of an element, this approach reduces to the simple idea of using two concentrated sources at either end of that edge [78].
• A third formulation is to separate the field in two parts – one part is a standard field produced by an ideal dipole in an infinite homogeneous domain and the other part is a solution in the closed sourceless domain under boundary conditions that correct the current movement across boundaries between regions of different conductivity [78].
• In the Laplace formulation, a small volume containing the dipole is removed and fixed boundary conditions are applied at all nodes on the surface of the removed volume. This can be interpreted as replacing current sources by an estimate of the equivalent voltage sources [78].
• A fifth formulation is the blurred dipole model, where source and sink monopoles are divided over the neighbouring nodes. In most cases the source and sink monopoles do not coincide with nodes of the FEMmesh. Therefore a way to represent the dipole is by a summation of monopoles placed at neighbouring nodes [79].
A comparison of the resulting surface potentials using the first four methods with the exact analytical solution using ideal dipoles (with an infinitesimal separation between the two poles, an infinite total current exiting one pole and entering the other, and a finite dipole moment, which is the product of the current and separation) in a 4layer concentric sphere was made in [78]. It was found that the third formulation gives the best performance for both transverse and radial dipoles (followed by the Laplace formulation for radial dipoles).
A recent innovation [80, 81] is to consider current monopoles (point sources/sinks) instead of dipoles. Using the equivalentcurrent inverse solution (ECS) approach for p grid locations, only p variables need to be determined in the inverse problem, whereas if a dipole is placed at each of the p grid locations, the solution space consists of 3p unknown variables because each dipole has 3 directional components. This results in an advantage of using current monopoles instead of current dipoles as demonstrated in [81] where it is shown that the time required to calculate the forward matrix in realistic finiteelement head models using the conventional approach may be reduced by one third. It is further shown in [80] that ECS imaging is equivalent to the equivalentdipole inverse solution (EDS) in that it provides sourcesink distribution corresponding to the dipole sources, but additional information is needed to determine the current flows (by combining ECS and EDS estimates).
In general the stiffness matrix is very big, making the computation of the electrode potentials very computationally intensive. To solve equation (42), iterative solvers for large sparse systems are used as given in [82]. Some techniques have been proposed to reduce the computational burden and increase efficiency as will be illustrated in section 5. A freely licensed software package that implements both FEM end BEM is NEUROFEM [79, 83–85].
The finite difference method
Isotropic media (iFDM)
For volumes G, which contain a current monopole, I_{ P }becomes I or I. α_{ i }has the dimension of Ω^{1} and corresponds with the conductance between P and Q_{ i }. Furthermore, for I_{ P }= 0 Kirchoff's law holds at the node P. For each node of a cubic grid we obtain a linear equation given by (44). The unknown potentials at the n computational points are represented by V ∈ ℝ^{n×1}. The source terms represented by I ∈ ℝ^{n×n}are calculated in each of the n cubes utilizing equation (45). Notice that in the linear equation (44) only the neighbouring computational points are included. The system matrix A ∈ ℝ^{n×n}has at most six offdiagonal elements and is a sparse matrix. In matrix notation one can write:
A·V = I
To solve this large sparse set of equations iterative methods are used. A discussion of the most popular solvers will be discussed in section. More extensive literature on this method can be found in [29, 87–92].
The finite difference method in anisotropic media (aFDM)
The differential equation (7) with Dirichlet and Neumann boundary conditions can be transformed into a set of linear equations even in the case of anisotropic media. This approach uses a cubic grid in which each cube (or element) has a conductivity tensor. In anisotropic tissues the conductivity tensor can vary between neighbouring elements. There are two ways to have anisotropy. In general, the directions of the anisotropy can be in any direction. Using tensor transformations, the matrix representation of the concuctivity tensor can be deduced. In a more specific case, the directions of the anisotropy are limited along the axes of the coordinate system of the headmodel. In this case, the anisotropy is orthotrophic [93]. In the next paragraphs, the general case as shown in [94] is depicted.
where ${\sigma}_{i}^{j}$ are the conductivities in the principal directions at element j, respectively. The matrix representation has to be transformed to a global cartesian coordinate system of the head. Therefore a rotation matrix has to be applied. The matrix representation of the conductivity tensor at element j in the cartesian system of the head is then given by ${\sigma}_{head}^{j}={T}_{j}^{T}{\sigma}^{j}{T}_{j}$, where T is a rotation transfer matrix from the local coordinate system to the global coordinate system [95] and ^{ T }denotes the transpose of a matrix.
A finite difference method which can handle anisotropic properties of tissues was presented by Saleheen et al. [96]. Here, the authors used a transition layer technique [97] to derive a finite difference formulation for the Laplace's equation that is valid everywhere in a piecewise inhomogeneous anisotropic medium, where the conductivity tensor can have a different value in each element of the grid. This formulation is extended to Poisson's equation by Hallez et al. [94].
For nodes at the corners of the compartments as illustrated in figure 12, for example node 11, the boundary normal cannot be obviously defined. Therefore, the Neumann boundary equations (12) and (11) contain singularities in spatial derivatives of the conductivities. The method presented in [94] and [96] has an advantage if one wants to enforce such a Neumann boundary condition: the formulation allows a discrete change or discontinuity in conductivity between neighbouring elements and will automatically incorporate the boundary between two different materials. In short, the boundary condition is already implicitly formulated in equation (47).
For each node a linear equation can be written as in equation (47), and for all computational points a set of linear equations is obtained: A·V = I. Due to the large size of the system, iterative methods have to be used.
Comparing the various numerical methods
A comparison of the different methods for solving Poisson's equation in a realistic head model is presented.
BEM  FEM  iFDM  aFDM  

Position of computational points  surface  volume  volume  volume 
Free choice of computational points  yes  yes  no  no 
System matrix  full  sparse  sparse  sparse 
Solvers  direct  iterative  iterative  iterative 
Number of compartments  small  large  large  large 
Requires tesselation  yes  yes  no  no 
Handles anisotropy  no  yes  no  yes 
A first difference between BEM and FEM or FDM is the domain in which the solutions are calculated. In the BEM the solutions are calculated on the boundaries between the homogeneous isotropic compartments while in the FEM and FDM the solution of the forward problem is calculated in the entire volume. Subsequently, the FEM and FDM lead to a larger number of computational points than the BEM. On the other hand, the potential at an arbitrary point can be determined with FEM and FDM by interpolation of computational points in its vicinity, while for the BEM it is necessary to reapply the Barnard formula [62] and numerical integration.
Another important aspect is the computational efficiency. In the BEM, a full matrix (I  C), represented in equation (31), needs to be inverted. When the scalp potentials need to be known for another dipole, V_{ 0 }in equation (31) needs to be recalculated and multiplied with the already available (I  C)^{1}. Hence once the matrix is inverted, only a matrix multiplication is needed to obtain the scalp potentials. This limited computational load is an attractive feature when solving the inverse problem, where a large number of forward evaluations need to be performed.
For the FEM and the FDM, a direct inversion of the large sparse matrices found in (42) and (46) is not possible due to the dimension of the matrices. Typically at least 500,000 computational points are considered which leads to system matrices of 500,000 equations with 500,000 unknowns which cannot be solved in a direct manner with the computers now available. However matrices found in FEM and FDM can be inverted for a given source configuration or righthand side term, utilizing iterative solvers such as the successive overrelaxation method, the conjugate gradient method [82], or algebraic multigrid methods [98, 99] (see next section). A disadvantage of the iterative solvers is that for each source configuration the solver has to be reapplied. The FEM and FDM would be computationally inefficient when for each dipole an iterative solver would need to be used. To overcome this inefficiency the reciprocity theorem is used, as will be explained in section.
When a large number of conducting compartments are introduced, a large number of boundaries need to be sampled for the BEM. This leads to a large full system matrix, thus a lower numerical efficiency. In FEM and FDM modeling, the heterogeneous nature of realistic head models will make the stiffness matrix less sparse and badly conditioned. Moreover, the incorporation of anisotropic conductivities will decrease the sparseness of the stiffness matrix. This can lead to an unstable system or very slow convergence if iterative methods are used. To obtain a fast convergence or a stable system, preconditioning should be used. Preconditioning transforms the system of equations Ax = b into a preconditioned system M^{1}Ax = M^{1}b, which has the same solution as the orignal system. M is a preconditioning matrix or a preconditioner and its goal is to reduce the condition number (ratio of the largest eigenvalue to the smallest eigenvalue) of the stiffness matrix towards the optimal value 1. Basic preconditioning can be used in the form of Jacobi, GaussSeidel, Successive OverRelaxation (SOR) and Symmetric Successive OverRelaxation (SSOR). These are easily implemented [100]. More advanced methods use incomplete LU factorization and polynomial preconditioning [93, 100].
For the FDM in contrast with the BEM and FEM, the computational points lie fixed in the cube centers for the isotropic approach and at the cube corners for the anisotropic approach. In the FEM and BEM, the computational points, the vertices of the tetrahedrons and triangles, respectively, can be chosen more freely. Therefore, the FEM can better represent the irregular interfaces between the different compartments than the FDM, for the same amount of nodes. However, the segmented medical images used to obtain the realistic volume conductor model are constructed out of cubic voxels. It is straightforward to generate a structured grid used in FDM from these segmented images. In the FEM and the BEM, additional tessellation algorithms [101] need to be used to obtain the tetrahedron elements and the surface triangles, respectively.
Finally, it is known that the conductivities of some tissues in the human head are anisotropic such as the skull and the white matter tissue. Anisotropy can be introduced in the FEM [102] and in the FDM [96], but not in the BEM.
Reciprocal approaches
In the literature one finds two approaches to solve the forward problem. In the conventional approach, the transfercoefficients making up the matrix G in equation (17) are obtained by calculating the surface potentials from dipole sources via Poisson's equation. The calculations are made for each dipole position within the head model and the potentials at the electrode positions are recorded.
In the reciprocal approach [103], Helmholtz' principle of reciprocity is used. The electric field that results at the dipole location within the brain due to current injection and withdrawal at the surface electrode sites is first calculated. The forward transfercoefficients are obtained from the scalar product of this electric field and the dipole moment. Calculations are thus performed for each electrode position rather than for each dipole position. This speeds up the time necessary to do the forward calculations since the number of electrodes is much smaller than the number of dipoles.
The general idea of reciprocity
Mathematical treatment
A mathematical treatment for a digitized volume conductor model is considered. Consider a digitized volume conductor model with n computational points or nodes. At each of the nodes the potential V_{ i }with i = 1...n is calculated for given sources which are the current monopoles I_{ i }with i = 1...n. Poisson's equation can then be transformed to a linear equation at each node, as illustrated for the FEM and FDM in subsections and. This set of linear equations can be written in matrix notation. The system matrix then becomes A ∈ ℝ^{n×n}and has the following properties: it is sparse, symmetric and regular. One can write:
A·V = I,
with A_{*∘} the minor for row * and column ∘.
Furthermore, A_{*∘} is equal to A_{∘*} due to the fact that A is symmetric. Hence, (eqn.(49)  eqn.(50))/I_{ f }equals (eqn.(51)  eqn.(52))/I_{ k }. Subsequently the reciprocity theorem is deduced:
I_{ k }(V_{ k } V_{ l }) = I_{ f }(V_{ f } V_{ g }).
Reciprocity for a dipole source with random orientation
In a similar way, U_{ AB }can be calculated for a dipole located at r oriented along the yaxis and the zaxis.
with ∇V(r) = (∂V(r)/∂x, ∂V(r)/∂y, ∂V(r)/∂z)^{ T }∈ ℝ^{3×1}.

A fictive current I_{ AB }of arbitrary value is introduced which enters the head at electrode A and leaves the head at electrode B.

Utilizing the FDM the potentials V(hi, hj, hk) can be calculated with h the internode spacing and i, j, k the node numbers along the Cartesian axes. Figure 15 illustrates the equipotential lines and current density vectors J = σ∇V in the brain region, with ∇V = (∂V/∂x, ∂V/∂y, ∂V/∂z)^{ T }. The partial derivative ∂V/∂x is approximated by [V(h(i + 1), hj, hk)  V(h(i  1), hj, hk)]/2h. The partial derivatives ∂V/∂y, ∂V/∂z are obtained in a similar way.

U_{ AB }the potential difference between the scalp electrodes A and B generated by the dipole at position r and dipole moment d is obtained by applying eqn. (55). When r does not coincide with a node, then ∇V(r) is obtained with trilinear interpolation [104].
By solving only one forward calculation numerically, by introducing current monopoles at electrodes A and B, and storing the obtained node potentials in a data structure, U_{ AB }is obtained for every dipole position and orientation.
If N scalp electrodes are used to measure the EEG, N  1 electrode pairs can be found with linear independent potential differences. Therefore N  1 numerical forward calculations are performed and stored in data structures. In addition, the N  1 potential differences at the N  1 electrode pairs are transformed in N average referenced potentials at the N electrodes.
Reciprocity has been applied in the literature in conjunction with BEM [105], FEM [106] and FDM [29, 89, 91, 94, 107].
Solving large sparse linear systems applied in FEM and FDM
Properties of the system matrix
If the linear system resulting is rewritten from equations (41), (44) and (47) in algebraic form as Ax = b, the system matrix A = {a_{ ij }} has the following properties.
From the coefficients in the linear equations one can see that the coefficient connecting a computational point V_{ l }to a neighbouring point V_{ k }is identical to the coefficient connecting V_{ k }to V_{ l }, thus A is symmetric. Moreover for FEM, it can be shown that the stiffness matrix A is a symmetric positive definite matrix [100].
From equations (44) and (47) it is easy to see that the righthand side of our problem satisfies this condition. The vector, b, represent a dipole or a multipole, hence the sum of the elements are zero. Thus the problem Ax = b possesses infinitely many solutions differing only in an additive constant.
Instead of solving the singular linear system, another possibility is to transform it into a regular one and solve this instead. The regular problem is chosen such that its unique solution belongs to the set of solutions of the original singular system. The easiest approach is to fix the value of a computational point to 0. The special structure of the matrix then allows us to cancel the corresponding row and column in A and also the respective entry in the righthand side vector b. This leads to a problem with a regular matrix and its solution obviously solves the initial problem with the particular computational point set zero.
Another important aspect of the matrix A is its sparseness. Every matrix row contains a few nonzero offdiagonal elements. This leads to a very small ratio of nonzero to overall entries resulting in a very sparse matrix.
Iterative solvers
The following methods will be discussed:

Successive overrelaxation (SOR)

Conjugate gradients (CG)

Preconditioned conjugate gradient method (PCG)

Algebraic multigrid (AMG)
While these methods have been developed for regular linear systems, they can also be applied in our semidefinite case. In the case of a consistent righthand side, semiconvergence can be guaranteed for SOR and (P)CG, while for AMG theoretical results are more complicated [108]. A summary of each method is given based on [104] for the first three and [109, 110] for the last method.
Successive overrelaxation
The SOR method is a representative of classical stationary methods. It is known to be a nonoptimal choice as far as convergence is concerned, but has a very simple structure. Thus it is a good candidate for an optimised implementation.
A linear system of equations Ax = b,
a_{i1}x_{1} + ⋯ + a_{ ii }x_{ i }+ ⋯ + a_{ in }x_{ n }= b_{ i }, i = 1,...,n,
where ρ(B) is the spectral radius or the maximum of the absolute eigenvalues of the Jacobi iteration matrix. During the SOR procedure, the ω can be altered using this formula to obtain a faster convergence [100].
Conjugate gradients
The CG method is the typical algorithm and is especially suited for symmetric positive definite matrices, for which it was originally devised. CG is a descendant of the method of steepest descent, that avoids repeated search in the same directions by making search directions orthogonal to each other in the energy (L2) norm associated with the matrix.
In the CG method the iterate x^{(k+1) }is computed via
x^{(k+1) }= x^{(k) }+ α^{(k)}d^{(k)},
The first search direction is just the residual of the initial guess d^{(0)} = r^{(0)}, where the residual is defined by r^{(k) }= b  Ax^{(k)}. The kth search direction is computed from the previous one via
d^{(k) }= r^{(k) }+ β^{(k)}d^{(k1)},
Preconditioned conjugate gradients
The convergence of the CG method depends on the condition of the problem matrix. More precisely it is the distribution of the eigenvalues of the matrix that determines the convergence. The distribution of the eigenvalues is also known as the spectrum of a matrix. Loosely speaking, the more eigenvalues lie close together in clusters, the faster the convergence. If the eigenvalues are widely scattered, a situation one can often find for problems with a large number of unknowns, the convergence will be slow. CG is therefore seldom used without preconditioning. By preconditioning, the spectral properties of the linear system should be improved. Instead of solving Ax = b, one solves the modified system $\tilde{A}\tilde{x}=\tilde{b}$ with $\tilde{A}$ = E^{1}A(E^{1})^{ T }, $\tilde{x}$ = E^{ T }x and $\tilde{b}$ = E^{1}b. In practice, the matrix $\tilde{A}$ is never explicitly formed. Instead in each step of the PCG algorithm a linear system of the form Mz = r must be solved with the matrix M = EE^{ T }and with r the residual. Compare the pseudocode of the preconditioned CG in figure 17.
Algebraic multigrid
Comparing the iterative solvers
For the iFDM, the following conclusions by [111] were drawn from comparing the iterative solvers: (a) the algebraic multigridbased solver performed the task 1.8–3.5 times faster, platform depending, than the secondbest contender, (b) there is no need to introduce a reference potential which forces a unique solution and (c) neither the grid nor matrixbased implementation of the solvers consistently gives a smaller run time.
Wolters et al. [79] investigated the parallel implementation of iterative solvers. If the JacobiCG (Jacobi Preconditioned Conjugate Gradient) solver on a single processor is taken as a reference, a speedup of 75 for a realistically shaped high resolution head model is achieved with the parallel AMGCG (Algebraic multiGrid preconditioned congjugate gradient) solver on 12 processors. This speedup can be attributed to the algebraic multigrid preconditioning (speedup of 7.5), and to the parallelization on 12 processors (speedup of 10). The required relative solution accuracy was 10^{8}.
Summary
The aim of this work was to give newcomers in the field of EEG source localization an overview of the stateoftheart methods applied to solve the forward problem. Multiple references to the work of authors active in this area were provided.
The postsynaptic potentials at apical dendrites of pyramidal cells are suggested to be the generators of the EEG. The extracellular electric currents generated by these cells obey the quasistatic conditions, i.e. all currents behave as if they were stationary at each instance in time. The electrical conductivity of a tissue can be isotropic, identical in all directions (e.g. fat, cerebrospinal fluid), or, anisotropic, not identical in all directions (e.g. white matter and skull). For both cases Poisson's equation with the Neumann and Dirichlet boundary conditions was derived. The active cluster of pyramidal cells was modelled with a current dipole.
Finding the electrode potentials for a given dipole source configuration is solving the socalled forward problem. The first models used were threeshell spherical head models. Analytical solutions exist here to solve the forward problem. To have a more accurate resolution, realistically shaped head models need to be constructed. These models can be obtained by segmenting MR/CT images to extract different conducting compartments associated with certain tissues. White matter anisotropy can be obtained from MR diffusion tensor images.
Various numerical methods can be used to solve the forward problem in a realistically shaped head model. BEM, FEM, iFDM and aFDM were discussed. For BEM the computational points are located on the surfaces between isotropic conducting compartments while for the other methods the computational points are located in the entire volume. Furthermore for BEM and FEM the computational points can be chosen freely. One could, for example, place more points in areas where more irregular shapes occur. An additional tessellation algorithm to position the computational points is then required. For FDM the cubic grid is rigid. This gives the user the opportunity to import directly from 3D medical images where cubical voxels are also used. Note also that for both FEM and aFDM anisotropic compartments can be used.
The reciprocity theorem was introduced to speed up the time necessary to solve the forward problem. The electric field that results at the dipole location within the brain due to current injection and withdrawal at the surface electrode sites is first calculated. The forward transfercoefficients are obtained from the scalar product of this electric field and the dipole moment. Calculations are thus performed for each electrode position rather than for each dipole position. This speeds up the time necessary to perform the forward calculations since the number of electrodes is much smaller than the number of dipoles.
For FEM and FDM a large linear system is generated with a sparse system matrix and a righthand side representing the electrical sources. Solving the forward problem is solving this linear system. Direct solvers cannot be used as the number of unknowns, the potentials at the computational points, is too large. Here iterative solvers for large sparse linear systems were utilized. The unknowns are only calculated for a given righthand side. The successive overrelaxation method, the conjugate gradient method, the preconditioned conjugate gradient method, and the multigrid method were discussed. The last method was the most promising computationtime wise.
Discussion and new trends
In this section interesting/necessary evolutions with respect to the forward problem in EEG source localization will be discussed. The following topics are raised: (a) a promising way to obtain tissue conductivity by magnetic resonance electromagnetic impedance tomography (MREIT); (b) combined EEG/MEG source localization; (c) incorporating invasive electrodes to overcome the disadvantages of the skull; (d) dipole localization benefits of improving the SNR of the EEG by blind source separation techniques; (e) combining EEG with functional magnetic resonance imaging (fMRI) to yield more accurate localization; (f) the necessity for a grand benchmark study to compare the performance of the different numerical methods on the same dataset; (g) use of advanced numerical approaches of the FEM and/or FDM; (h) numerical approaches for dipole modelling in the forward problem.
(a) One of the main problems in EEG source localization is the uncertainty of the tissue conductivity. Although there are a lot of studies concerning this topic, the actual conductivity is not well established (and may change from person to person with age). In section 3, a distinct set of brain tissues and assigned bulkconductivity to each of them was explained, thus obtaining a piecewise homogeneous head model. In reality the conductivity within brain tissue is place dependent and, thus, variable. The boundaries of brain tissues are in reality not discrete but a continuum. The technique to measure conductivity at a specific place in the brain is impedance tomography. A recent promising technique, MREIT, utilizes the information in both magnetic as well as electric fields (induced by injection current at the surface electrodes) to build a conductivity profile of the human head in 3D. While EIT is limited by the boundary measurements of currentvoltage data, MREIT utilizes the internal magnetic flux density data obtained using a Magnetic Resonance Imaging (MRI) scanner. When a current I is injected into an electrically conducting body through a pair of surface electrodes, an electric current density J and a magnetic field density B are formed insed the conducting body. The magnetic field density B can be measured inside a MR scanner, J can be calclulated by J = ∇ × B/μ_{0}, μ_{0} is the magnetic permeability in the free space. The conductivity image reconstruction problem in MREIT means to find a conductivity map σ from the data set [I, J_{ m }, V_{ m }] → σ, where J_{ m }is the measured current density inside the subject, V_{ m }is the measured voltage between the electrodes, and σ is the conductivity image to be reconstructed [112]. However, the technique is highly sensitive to noise and there is an open problem on the uniqueness of the reconstructed conductivity image. To incorporate anisotropic conductivity reconstruction one should use in addition diffusion tensor images [113]. This way the conductivity can be measured very locally [114, 115]. This could not only help us establish an accurate volume conductor model but could also give us the reciprocal current sources by measurment rather than numerical calculation.
(b) Magnetoencephalography or MEG measures the induction outside the head, generated by neuronal activity. Unlike EEG, MEG is less sensitive for the conductivity of the skull. However, MEG are unable to measure radially oriented dipoles. The MEG equipment consists of superconducting quantum interference devices (SQUIDs) that can measure very low variations of magnetic field differences. A major drawback of MEG is the huge sensitivity to instrumental and environment noise. The use of superconduting elements can minimize the instrumental noise. Noise from highfrequent electromagnetic waves, like radiowaves, can be eleminated using a shielded room. Slow electromagnetic waves, like passing cars, can be minimized by the use of gradiometers. These effort to minimize the noise sources, limit the possibility to conduct a longterm monitoring. Nevertheless, the combined use of EEG and MEG has shown to be beneficial for stimuluslocked brain activations. More recently, EEG and MEG have come to be viewed as complementary rather than competing modalities. The combine EEG/MEG measurement can compensate each one's limitations and can be a very succesful modality. An accurate modelling of the human head can improve the solution of both EEG and MEG [116, 117].
(c) Due to the lower conductivity of the skull compared to the surrounding tissue, the potentials generated by a source in the brain are smeared out over the scalp surface, the skull is acting as a spatially lowpass filter. In some epilepsy patients, depth and grid electrodes are implanted in their brain. While grid electrodes are arranged in an array of 8, 16 or 32 electrodes and measure the electrical activity at the cortical surface of the brain, depth electrodes are implanted in the brain near the presumably active brain structures. These electrodes measure the electrical activity without the shielding effect of the skull. Initial studies have shown that including this information in source localization may improve the accuracy [118]. Although, the brain cannot fully be surrounded by grid electrodes, as the surgery of the patient would be to intensive. As the signaltonoise ratio of the grid electrodes is larger than scalp electrodes, it is difficult to use both grid electrodes and scalp electrodes at the same time in the dipole estimation problem. One way to circumvent this is to create an a priori distribution of the brain activity using the grid electrodes and then use the scalp electrodes to do a more precise estimation. In the same aspect, a multipole model as a source model for intracranial EEG can become important as two dipoles with opposite orientation travel along the axon very close to each other.
(d) Noise coming from EEG background activity (i.e. other brainwaves than the ones you are interested in), artefacts from extracerebral sources (such as eye movements and muscle activation) and instrument quantization noise are inherent to the EEG and limit the accuracy of source localization. Removing the noise from the EEG signal is important and should be investigated. It is not necessary to incorporate anisotropic compartments when knowing that the dipole error due to noise in the EEG is of a higher magnitude. New algorithms to filter noise and artefacts from the EEG are Blind Source Separation (BSS) techniques like Independent Component Analysis (ICA) and Canonical Correlation Analysis (CCA). For more detail on these issues we refer to the literature [119–121].
(e) While the spatial resolution of the EEG is low because of the noise of the EEG signals, other modalities have a high spatial resolution. Functional magnetic resonance imaging measures locally the blood oxygen level in the brain. The fMRI scan highlights the activity present at a specific brain area. This technique has a high spatial resolution, but (unlike the EEG) has a low temporal resolution. Using both modalities, EEG and fMRI, in source localization should improve the accuracy [116, 122–124].
(f) The field of forward modelling has grown extensively in the last years. Now more than ever there is need for benchmark studies to compare different numerical techniques. A conceptual benchmarking study has been done by Pruis et al. [125], but lacks a quantitative assessment of the pros and cons of the numerical methods. A benchmarking study can investigate which numerical method is preferable to others, in a specific situation. This needs an extensive and diverse database of head models accessible for researchers.
(g) Making FEM or FDM more accurate can be done in 2 ways: (i) increasing the number of nodes and (ii) increasing the order of an element. Increasing de number of nodes is equivalent to reducing element size and is know as classical (or the hversion) FEM, while increasing the order of an element is known as the pversion of FEM. The generalisation of both variants is the hpversion of FEM (hpFEM) and is a spectral method where the convergence is achieved by simultaneously refining the mesh and increasing the approximation order. It originated in the late 1980s and is becoming very important due to its potential of unconditional exponential convergence [126, 127]. The main idea is to place large high order elements in regions where the potential is smooth, and small elements with low order in regions close to singularities. One could obtain the accuracy of the system but decrease significantly the number of elements (degrees of freedom) and also the computation time [128]. However, the usefullness of hpFEM in modeling the forward problem in EEG source localization is yet to be established.
(h) In the forward problem, the mathematical dipole is a current source and sink infinitesimally close to each other. This introduces a singularity on the righthand side of the Poisson equation that has to be treated specifically. Several approaches exist to deal with this singularity. These approaches were summarized and benchmarked in [129]. In this article, it was shown that the potential distribution changes significantly when using different approaches for the dipole model. Further research in this field can improve our understanding in numerical methods to solve the forward problem and give insights towards correct numerical modelling of the dipole.
Declarations
Acknowledgements
Hans Hallez is funded by a grant of the Institute for the Promotion of Innovation through Science and Technology in Flanders (IWTVlaanderen). Bart Vanrumste is postdoctoral fellow at K. U. Leuven in Belgium and funded by the 'Programmatorische Federale Overheidsdienst Wetenschapsbeleid' of the Belgian Government. Research of Sabine Van Huffel, Anneleen Vergult, Wim De Clercq and Bart Vanrumste is supported by Research Council KUL: GOAAMBioRICS, CoE EF/05/006 Optimization in Engineering, IDO 05/010 EEGfMRI; Flemish Government: FWO: projects, G.0360.05 (EEG, Epileptic), FWOG. 0321.06 (Tensors/Spectral Analysis), G.0341.07 (Data fusion), research communities (ICCoS, ANMMM); IWT: PhD Grants; Belgian Federal Science Policy Office IUAP P6/04 (Dynamical systems, control and optimization, 2007–2011); EU: BIOPATTERN (FP62002IST 508803). Joseph Muscat, Roberta Grech, Simon G. Fabri and Kenneth P. Camilleri would like to acknowledge that part of this work was supported by University of Malta Grant LBA73695.
Authors’ Affiliations
References
 Berger H: Über das Elektroenkephalogramm des Menschen. Archiv für Psychiatrie und Nervenkrankheiten. 1929, 87: 527570. 10.1007/BF01797193.Google Scholar
 Niedermeyer E, Lopes Da Silva F: Electroencephalography. 1993, Baltimore: Williams and WilkinsGoogle Scholar
 Baillet S, Mosher JC, Leahy RM: Electromagnetic Brain Mapping. IEEE Signal Processing Magazine. 2001, November: 1430. 10.1109/79.962275.Google Scholar
 Kiloh LG, Mc Comas AJ, Osselton JW, Upton ARM: Clinical Electroencephalography. 1981, ButterworthsGoogle Scholar
 Gray H: Gray's Anatomy. 1973, Longman Group LtdGoogle Scholar
 Gulrajani RM: Bioelectricity and Biomagnetism. 1998, New York: John Wiley & Sons, IncGoogle Scholar
 Johnston D, Wu SMS: Foundations of Cellular Neurophysiology. 1994, the MIT Press, Massachusetts Institute of TechnologyGoogle Scholar
 Gulrajani RM: Bioelectricity and Biomagnetism, chap Electroencephalography. 1998, John Wiley & Sons, Inc, 469524.Google Scholar
 Schaul N: The Fundamental Neural Mechanisms of Electroencephalography. Electroencephalography and Clinical Neurophysiology. 1998, 106: 101107. 10.1016/S00134694(97)001119.PubMedGoogle Scholar
 Malmivuo J, Plonsey R: Bioelectromagnetism: Principles and Applications of Bioelectric and Biomagnetic Fields. 1995, New York: Oxford University PressGoogle Scholar
 Speckmann EJ, Elger CE: Introduction to the Neurophysiology, basis of the EEG and DC potentials. Electroencephalography, Basic Principles, Clinical Applications and Related Fields. Edited by: Niedermeyer E, Lopes da Silva F. 1987, Urban and Schwarzenberg, 113. 2Google Scholar
 Lopes da Silva F, Van Rotterdam A: Electroencephalography: Basic Principles, Clinical Applications and Related Fields, chap Biophysical aspects of EEG and magnetoencephalogram generation. 1993, Williams and Wilkins, 7891.Google Scholar
 Plonsey R, Heppner DB: Considerations of quasistationarity in electrophysiological systems. Bulletin of Mathematical Biophysics. 1967, 29 (4): 657664. 10.1007/BF02476917.PubMedGoogle Scholar
 Vanrumste B: EEG dipole source analysis in a realistic head model. PhD thesis. 2001, Ghent UniversityGoogle Scholar
 Plonsey R: Bioelectric Phenomena. 1969, McGrawHill, New YorkGoogle Scholar
 Rush S, Driscoll DA: Current distribution in the brain from surface electrodes. Anesth. Analgesia. 1968, 47: 717723.Google Scholar
 Wolters C: Influence of Tissue Conductivity Inhomogeneity and Anisotropy on EEG/MEG based Source Localisation in the Human Brain. PhD thesis. 2003, Universität LeipzigGoogle Scholar
 Wolters CH, Anwander A, Tricoche X, Weinstein D, Koch MA, MacLeod RS: Influence of tissue conductivity anisotropy on EEG/MEG field and return current computation in a realistic head model: a simulation and visualization study using highresolution finite element modeling. Neuroimage. 2006, 30 (3): 813826. 10.1016/j.neuroimage.2005.10.014.PubMedGoogle Scholar
 Nicholson P: Specific impedance of cerebral white matter. Experimental Neurology. 1965, 13: 386401. 10.1016/00144886(65)901263.PubMedGoogle Scholar
 Geddes LA, Baker LE: The specific resistance of biological material – a compendium of data for the biomedical engineer and physiologist. Med Biol Eng. 1967, 5 (3): 271293. 10.1007/BF02474537.PubMedGoogle Scholar
 Basser P, Mattiello J, LeBihan D: MR Diffusion Tensor Spectroscopy and Imaging. Biophysical Journal. 1994, 66: 259267.PubMed CentralPubMedGoogle Scholar
 Tuch D, Wedeen V, Dale A, George J, Belliveau J: Conductivity tensor mapping of the human brain using Diffusion Tensor MRI. Proceedings in National Academie of Science. 2001, 98 (20): 1169711701. 10.1073/pnas.171473898.Google Scholar
 Haueisen J, Tuch D, Ramon C, Schimpf P, Wedeen V, George J, Belliveau J: The Influence of Brain Tissue Anisotropy on Human EEG and MEG. NeuroImage. 2002, 15: 159166. 10.1006/nimg.2001.0962.PubMedGoogle Scholar
 Hallez H, Van Hese P, Vanrumste B, Boon P, D'Asseler Y, Lemahieu I, Van de Walle R: Dipole Localization Errors due to not Incorporating Compartments with Anisotropic Conductivities: Simulation Study in a Spherical Head Model. International Journal of Bioelectromagnetism. 2005, 7: 134137.Google Scholar
 Pursula A, Nenonen J, Somersalo E, Ilmoniemi E, Katila T: Bioelectromagnetic calculations in anisotropic volume conducters. Proceedings of Biomag2000. 2000, 659662.Google Scholar
 Baumann S, Wozny D, Kelly S, Meno F: the Electrical Conductivity of Human Cerebrospinal Fluid at Body Temperature. IEEE Transactions on Biomedical Engineering. 1997, 44 (3): 220223. 10.1109/10.554770.PubMedGoogle Scholar
 Awada K, Jackson D, Baumann S, Williams J, Wilton D, Fink PW, Prasky BR: Effect of conductivity uncertainties and modeling errors on EEG source localization using a 2Dmodel. IEEE Transactions on Biomedical Engineering. 1998, 45 (9): 11351145. 10.1109/10.709557.PubMedGoogle Scholar
 Oostendorp T, Delbeke J, Stegeman D: The conductivity of the human skull: results of in vivo and in vitro measurements. Biomedical Engineering, IEEE Transactions on. 2000, 47 (11): 14871492. 10.1109/TBME.2000.880100.Google Scholar
 Vanrumste B, Van Hoey G, Van de Walle R, D'Havé M, Lemahieu I, Boon P: Dipole location errors in electroencephalogram source analysis due to volume conductor model errors. Medical and Biological Engineering and Computing. 2000, 38 (5): 528534. 10.1007/BF02345748.PubMedGoogle Scholar
 Gonçalves S, de Munck J, Verbunt J, Bijma F, Heethaar R, Lopes da Silva F: In vivo measurement of the brain and skull resistivities using an EITbased method and realistic models for the head. Biomedical Engineering, IEEE Transactions on. 2003, 50 (6): 754767. 10.1109/TBME.2003.812164.Google Scholar
 Clerc M, Badier JM, Adde G, Kybic J, Papadopoulo T: Boundary element formulation for electric impedance tomography. ESAIM: Proceedings. 2005, 14: 6371.Google Scholar
 Gutirrez D, Nehorai A, Muravchik CH: Estimating brain conductivities and dipole source signals with EEG arrays. IEEE Trans Biomed Eng. 2004, 51 (12): 21132122. 10.1109/TBME.2004.836507.Google Scholar
 Lai Y, van Drongelen W, Ding L, Hecox KE, Towle VL, Frim DM, He B: Estimation of in vivo human braintoskull conductivity ratio from simultaneous extra and intracranial electrical potential recordings. Clin Neurophysiol. 2005, 116 (2): 456465. 10.1016/j.clinph.2004.08.017.PubMedGoogle Scholar
 Hoekema R, Wieneke GH, Leijten FSS, van Veelen CWM, van Rijen PC, Huiskamp GJM, Ansems J, van Huffelen AC: Measurement of the conductivity of skull, temporarily removed during epilepsy surgery. Brain Topogr. 2003, 16: 2938. 10.1023/A:1025606415858.PubMedGoogle Scholar
 Rojansky : Electromagnetic fields and waves. 1971, DoverGoogle Scholar
 de Munck J, Van Dijk B, Spekreijse H: Mathematical Dipoles are Adequate to Describe Realistic Generators of Human Brain Activity. IEEE Transactions on Biomedical Engineering. 1988, 35 (11): 960965. 10.1109/10.8677.PubMedGoogle Scholar
 He B, Yao D, Lian J: Highresolution EEG: on the cortical equivalent dipole layer imaging. Clinical Neurophysiology. 2002, 113 (2): 22735. 10.1016/S13882457(01)007404.PubMedGoogle Scholar
 Hara J, Musha T, Shankle WR: Approximating Dipoles from Human EEG Activity: The Effect of Dipole Source Configuration on Dipolarity Using Single Dipole Models. IEEE Transactions on Biomedical Engineering. 1999, 46 (2): 125129. 10.1109/10.740874.PubMedGoogle Scholar
 Karp PJ, Katila TE, Saarinen M, Siltanen P, Varpula TT: The normal human magnetocardiogram. II. A multipole analysis. Circ Res. 1980, 47: 117130.PubMedGoogle Scholar
 Faugeras O, Clement F, Deriche R, Keriven R, Papadopoulo T, Roberts J, Vieville T, Devernay F, Gomes J, Hermosillo G, Kornprobst P, Lingrand D: The Inverse EEG and MEG Problems: The Adjoint State Approach I: The Continuous Case. Tech rep. 1999, INRIAGoogle Scholar
 Jerbi K, Baillet S, Mosher JC, Nolte G, Garnero L, Leahy RM: Localization of realistic cortical activity in MEG using current multipoles. Neuroimage. 2004, 22 (2): 779793. 10.1016/j.neuroimage.2004.02.010.PubMedGoogle Scholar
 Jerbi K, Mosher JC, Baillet S, Leahy RM: On MEG forward modelling using multipolar expansions. Physics in Medicine and Biology. 2002, 47 (4): 523555. 10.1088/00319155/47/4/301.PubMedGoogle Scholar
 Frank E: Electric Potential Produced by Two Point Current Sources in a Homogeneous Conduction Sphere. Journal of Applied Physics. 1952, 23 (11): 12251228. 10.1063/1.1702037.Google Scholar
 Ary JP, Klein SA, Fender DH: Location of Sources of Evoked Scalp Potentials: Corrections for Skull and Scalp Thicknesses. IEEE Transactions on Biomedical Engineering. 1981, BME28 (6): 447452. 10.1109/TBME.1981.324817.Google Scholar
 Salu Y, Cohen LG, Rose D, Sato S, Kufta C, Hallett M: An Improved Method for Localizing Electric Brain Dipoles. IEEE Transactions on Biomedical Engineering. 1990, 37 (7): 699705. 10.1109/10.55680.PubMedGoogle Scholar
 Vanrumste B, Van Hoey G, Van de Walle R, D'Havé M, Lemahieu I, Boon P: Comparison of performance of spherical and realistic head models in dipole localization from noisy EEG. Medical Engineering & Physics. 2002, 24: 403418. 10.1016/S13504533(02)00036X.Google Scholar
 Huerta MA, Gonzalez G: The Surface Potentials Produced by Electric Sources in Stratified Spherical and Prolate Spheroidal Volume Conductors. International Journal of Electronics. 1983, 54 (5): 657671. 10.1080/00207218308938765.Google Scholar
 de Munck J, Peters M: A Fast Method to Compute the Potential in the Multiphere Model. IEEE Transactions on Biomedical Engineering. 1993, 40 (11): 11661174. 10.1109/10.245635.PubMedGoogle Scholar
 Zhang Z: A fast method to compute surface potentials generated by dipoles within multilayer anisotropic spheres. Physics in Medicine and Biology. 1995, 40: 335349. 10.1088/00319155/40/3/001.PubMedGoogle Scholar
 Meijs P, Peters MJ: The EEG and MEG using a model of eccentric spheres to describe the human head. IEEE Transactions on Biomedical Engineering. 1987, 34: 913920. 10.1109/TBME.1987.325929.PubMedGoogle Scholar
 Kariotou F: Electroencephalography in ellipsoidal geometry. Journal of Mathematical Analysis and Applications. 2004, 290: 32442. 10.1016/j.jmaa.2003.09.066.Google Scholar
 Vatta F, Bruno P, Inchingolo P: Multiregion BicentricSpheres Models of the Head for the Simulation of Bioelectric Phenomena. IEEE Transactions in Biomedical Engineering. 2005, 52 (3): 384389. 10.1109/TBME.2004.843258.Google Scholar
 Berg P, Scherg M: A Fast Method for Forward Computation of MultipleShell Spherical Head Models. Electroencephalography and Clinical Neurophysiology. 1994, 90: 5864. 10.1016/00134694(94)901139.PubMedGoogle Scholar
 Roth B, Gorbach A, Sato S: How well does a threeshell model predict positions of dipoles in a realistically shaped head?. Electroencephalography and Clinical Neurophysiology. 1993, 87: 175184. 10.1016/00134694(93)90017P.PubMedGoogle Scholar
 Roth BJ, Ko D, von AlbertiniCarletti IR, Scaffidi D, Sato S: Dipole Localization in Patients with Epilepsy Using the Realistic Shaped Head Model. Electroencephalography and Clinical Neurophysiology. 1997, 102 (3): 160166. 10.1016/S00134694(96)951115.Google Scholar
 Huiskamp G, Vroeijenstijn M, van Dijk R, Wieneke G, van Huffelen A: The Need for Correct Realistic Geometry in the Inverse EEG Problem. IEEE Transactions on Biomedical Engineering. 1999, 46 (11): 12811287. 10.1109/10.797987.PubMedGoogle Scholar
 Cuffin BN: Effects of Local Variations in Skull and Scalp Thickness on EEG's and MEG's. IEEE Transactions on Biomedical Engineering. 1993, 40: 4248. 10.1109/10.204770.PubMedGoogle Scholar
 Chauveau N, Franceries X, Doyon B, Rigaud B, Morucci JP, Celsis P: Effects of skull thickness, anisotropy, and inhomogeneity on forward EEG/ERP computations using a spherical threedimensional resistor mesh model. Human Brain Mapping. 2003, 21 (2): 8697. 10.1002/hbm.10152.Google Scholar
 Ermer JJ, Mosher JC, Baillet S, Leahy RM: Rapidly recomputable EEG forward models for realistic head shapes. Physics in Medicine and Biology. 2001, 46 (4): 12651281. 10.1088/00319155/46/4/324.PubMedGoogle Scholar
 He B, Musha T, Okamoto Y, Homma S, Nakajima Y, Sato T: Electric dipole tracing in the brain by means of the boundary element method and its accuracy. IEEE Transactions on Biomedical Engineering. 1987, 34 (6): 406414. 10.1109/TBME.1987.326056.PubMedGoogle Scholar
 Geselowitz DB: On Bioelectric Potentials in an Inhomogeneous Volume Conductor. Biophysics Jounal. 1967, 7: 111.Google Scholar
 Barnard ACL, Duck IM, Lynn MS: The Application of Electromagnetic Theory to Electrocardiography: I. Derivation of the Integral Equations. Biophysics Journal. 1967, 7 (5): 443462.Google Scholar
 Sarvas J: Basic Mathematical and Electromagnetic Concepts of the Biomagnetic Inverse Problem. Phys Med Biol. 1987, 32 (1): 1122. 10.1088/00319155/32/1/004.PubMedGoogle Scholar
 Barnard ACL, Duck IM, Lynn MS, Timlake WP: The Application of Electromagnetic Theory to Electrocardiography: II. Numerical Solution of the Integral Equations. Biophys J. 1967, 7 (5): 463491.PubMed CentralPubMedGoogle Scholar
 Meijs JWH, Weier OW, Peters MJ, van Oosterom A: On the Numerical Accuracy of the Boundary Element Method. IEEE Transactions on Biomedical Engineering. 1989, 36 (10): 10381049. 10.1109/10.40805.PubMedGoogle Scholar
 de Munck JC: A Linear Discretization of the Volume Conductor Boundary Integral Equation Using Analytically Integrated Elements. IEEE Transactions on Biomedical Engineering. 1992, 39 (9): 986990. 10.1109/10.256433.PubMedGoogle Scholar
 Frijns J, de Snoo S, Schoonhoven R: Improving teh accuracy of the boundary method by the use of secondorder interpolation functions. IEEE Transactions on Biomedical Engineering. 2000, 47: 13361346. 10.1109/10.871407.PubMedGoogle Scholar
 Gençer NG, Tanzer I: Forward problem solution of electromagnetic source imaging using a new BEM formulation with higherorder elements. Physics in Medicine and Biology. 1999, 44: 22752287. 10.1088/00319155/44/4/009.PubMedGoogle Scholar
 Lynn M, Timlake W: The use of multiple deflations in the numerical solution of singular systems of equations with application to potential theory. SIAM Journal of Numerical Analysis. 1968, 5 (2): 303322. 10.1137/0705027.Google Scholar
 Hämäinen M, Sarvas J: Realistic Conductivity Geometry Model of the Human Head for Interpretation of Neuromagnetic Data. IEEE Transactions on Biomedical Engineering. 1989, 36 (2): 165171. 10.1109/10.16463.Google Scholar
 Fuchs M, Drenckhahn R, Wischmann HA, Wagner M: An improved boundary element method for realistic volumeconductor modeling. IEEE Transactions on Biomedical Engineering. 1998, 45 (8): 980997. 10.1109/10.704867.PubMedGoogle Scholar
 Gençer N, AkalinAcar Z: Use of the isolated problem approach for multicompartment BEM models of electromagnetic source imaging. Physics in medicine and biology. 2005, 50: 30073022. 10.1088/00319155/50/13/003.PubMedGoogle Scholar
 AkalinAcar Z, Gençer N: An advanced boundary element method BEM implementation for the forward problem of electromagnetic source imaging. Physics in Medicine and Biology. 2004, 49: 50115028. 10.1088/00319155/49/21/012.PubMedGoogle Scholar
 Kybic J, Clerc M, Abboud T, Faugeras O, Keriven R, Papadopoulo T: A common formalism for the integral formulations of the forward EEG problem. IEEE Transactions on medical imaging. 2005, 24: 1228. 10.1109/TMI.2004.837363.PubMedGoogle Scholar
 Kybic J, Clerc M, Faugeras O, Keriven R, Papadopoulo T: Fast multipole acceleration of the MEG/EEG boundary element method. Physics in Medicine and Biology. 2005, 50: 46954710. 10.1088/00319155/50/19/018.PubMedGoogle Scholar
 Kybic J, Clerc M, Faugeras O, Keriven R, Papadopoulo T: Generalized head models for MEG/EEG: boundary element method beyond nested models. Physics in Medicine and Biology. 2006, 51: 13331346. 10.1088/00319155/51/5/021.PubMedGoogle Scholar
 Johnson CR: Numerical methods for bioelectric field problems. The biomedical engineering handbook. Edited by: Bronzino JD. 1995, CRC press, IEEE pressGoogle Scholar
 Schimpf P, Ramon C, Haueisen J: Dipole models for the EEG and MEG. IEEE Transactions on Biomedical Engineering. 2002, 49 (5): 409418. 10.1109/10.995679.PubMedGoogle Scholar
 Wolters C, Kuhn M, Anwander A, Reitzinger S: A parallel algebraic multigrid solver for finite element method based source localization in the human brain. Computing and Visualization in Science. 2002, 5: 165177. 10.1007/s0079100200980.Google Scholar
 He B, Yao D, Lian J, Wu D: An Equivalent Current Source Model and Laplacian Weighted Minimum Norm Current Estimates of Brain Electrical Activity. IEEE Transactions on Biomedical Engineering. 2002, 49 (4): 277288. 10.1109/10.991155.PubMedGoogle Scholar
 Muscat J, Grech R, Camilleri K, Fabri S, James C: An Improvement in EEG Forward Model Computational Efficiency. Proceedings of 2nd International Conference on Advances in Medical Signal and Information Processing. 2004, 99103.Google Scholar
 Datta BN: Numerical Linear Algebra and Applications. 1995, Brooks/Cole Publishing CompanyGoogle Scholar
 Wolters C, Grasedyck L, Hackbusch W: Efficient Computation of Lead Field Bases and Influence Matrix for the FEMbased EEG and MEG Inverse Problem. Inverse Problems. 2004, 20 (4): 10991116. 10.1088/02665611/20/4/007.Google Scholar
 Wolters C, Hartmann U, Koch M, Kruggel F: New Methods for Improved and Accelerated FEvolume conductor modelling in EEG/MEGSource Localisation. Proceedings of 4th Symposium on Computer methods in Biomechanics and Biomedical Engineering 99. 489494.Google Scholar
 IPNeurofem: A C++ class structured code for fast highresolution finite element method based EEG and MEG source analysis. [http://www.neurofem.com]
 Mitchell A, Griffiths D: The Finite Difference Method in Partial Differential Equations. 1980, John Willey and SonsGoogle Scholar
 Witwer JG, Trezek GJ, Jewett DL: The Effect of Media Inhomogeneities Upon Intracranial Electrical Fields. IEEE Transactions on Biomedical Engineering. 1972, BME19 (5): 352362. 10.1109/TBME.1972.324138.Google Scholar
 Marino F, Halgren E, Badier JM, Gee M, Nenev V: A Finite Difference Model of Electric Field Propagation in the Human Head: Implementation and Validation. Proceedings of the 19th Annual Northeast Bioengineering Conference. 1993, 8285.Google Scholar
 Laarne P, H E, Hyttinen J, Suihko V, Malmivuo J: Validation of a Detailed Computer Model for the Electric Fields in the Brain. Journal of Medical Engineering and Technology. 1995, 19 (2–3): 8487.PubMedGoogle Scholar
 Lemieux L, McBride A, Hand JW: Calculation of Electrical Potentials on the surface of a Realistic Head Model by Finite Differences. Physics in Medicine and Biology. 1996, 41: 10791091. 10.1088/00319155/41/7/001.PubMedGoogle Scholar
 Laarne P, Hyttinen J, Dodel S, Malmivuo J, Eskola H: Accuracy of two dipolar inverse algorithms applying reciprocity for forward calculation. Computers and Biomedical Research. 2000, 33 (3): 172185. 10.1006/cbmr.1999.1538.PubMedGoogle Scholar
 Laarne P, TenhunenEskelinen M, J H, Eskola H: Effect of EEG electrodes density on dipole localization accuracy using two realistically shaped skull resistivity models. Brain Topography. 2000, 12 (4): 249254. 10.1023/A:1023422504025.PubMedGoogle Scholar
 Neilson LA, Kovalyov M, Koles ZJ: A computationally efficient method for accurately solving the EEG forward problem in a finely discretized head model. Clin Neurophysiol. 2005, 116 (10): 23022314. 10.1016/j.clinph.2005.07.010.PubMedGoogle Scholar
 Hallez H, Vanrumste B, Van Hese P, D'Asseler Y, Lemahieu I, Van de Walle R: A finite difference method with reciprocity used to incorporate anisotropy in electroencephalogram dipole source localization. Physics in Medicine and Biology. 2005, 50: 37873806. 10.1088/00319155/50/16/009.PubMedGoogle Scholar
 Hinchey F: Vectors and Tensors for Engineers and Scientists. 1976, Wiley, John & Sons, Incorporated, [ISBN: 0470151943]Google Scholar
 Saleheen H, Kwong T: New Finite Difference Formulations for General Inhomogeneous Anisotropic Bioelectric Problems. IEEE Transactions on Biomedical Engineering. 1997, 44 (9): 800809. 10.1109/10.623049.PubMedGoogle Scholar
 Panizo M, Castellanos A, Rivas J: Finitedifference operators in inhomogeneouw anisotropic media. Journal of Applied Physics. 1977, 48 (3): 10541057. 10.1063/1.323779.Google Scholar
 Briggs WL: A multigrid tutorial. 1987, SIAMGoogle Scholar
 Hoekema R, Venner K, Struijk J, Holsheimer J: Multigrid solution of the potential field in modeling electrical nerve stimulation. Computers and Biomedical Research. 1998, 31: 348362. 10.1006/cbmr.1998.1486.PubMedGoogle Scholar
 Saad Y: Iterative methods for sparse linear systems. 2003, Philadelphia: SIAM, 2Google Scholar
 Thompson JF, Soni BK, Weatherrill NP: Handbook of grid generation. 1998, CRC PressGoogle Scholar
 Ottosen N, Peterson H: Introduction to the finite element method. 1992, Prentice hallGoogle Scholar
 Rush S, Driscoll DA: EEG Electrode Sensitivity – An Application of Reciprocity. IEEE Trans Biomed Eng. 1969, 16 (1): 1522.PubMedGoogle Scholar
 Press WH, Teukolsky SA, Vetterling WT, Flannery BP: Numerical Recipes in C. 1995, Cambridge University PressGoogle Scholar
 Finke S, Gulrajani R, Gotman J: Conventional and Reciprocal Approaches to the Inverse Dipole Localization Problem of Electroencephalography. IEEE Transactions on Biomedical Engineering. 2003, 50 (6): 657666. 10.1109/TBME.2003.812198.PubMedGoogle Scholar
 Weinstein D, Zhukov L, Johnson C: Leadfield bases for electroencephalography source imaging. Annals of biomedical engineering. 2000, 28: 10591065. 10.1114/1.1310220.PubMedGoogle Scholar
 Vanrumste B, Van Hoey G, Van de Walle R, D'Havé M, Lemahieu I, Boon P: The validation of the finite difference method and reciprocity for solving the inverse problem in EEG dipole source analysis. Brain Topography. 2001, 14 (2): 8392. 10.1023/A:1012909511833.PubMedGoogle Scholar
 Berman A, Plemmons R: Nonnegative Matrices in the Mathematical Sciences. No 9 in Classics in Applied Mathematics. 1994, SIAMGoogle Scholar
 Ruge JW, Stüben K: Algebraic multigrid (AMG). Multigrid Methods, Volume 3 of Frontiers in Applied Mathematics. Edited by: McCormick SF. 1987, SIAM, 73130.Google Scholar
 Briggs WL, Henson VE, McCormick SF: A Multigrid Tutorial. 2000, SIAM, 2Google Scholar
 Mohr M, Vanrumste B: Comparing iterative solvers for linear systems associated with the finite difference discretisation of the forward problem in electroencephalographic source analysis. Medical & Biological Engineering & Computing. 2003, 41: 7584. 10.1007/BF02343542.Google Scholar
 Oh SH, Han JY, Lee SY, Cho MH, Lee BI, Woo EJ: Electrical conductivity imaging by magnetic resonance electrical impedance tomography (MREIT). Magnetic Resonance in Medicine. 2003, 50: 875878. 10.1002/mrm.10588.PubMedGoogle Scholar
 Seo JK, Pyo PC, Chan Hyun, Kwon O, Woo EJ: Image reconstruction of anisotropic conductivity tensor distribution in MREIT: computer simulation study. Physics in Medicine and biology. 2004, 49: 43714382. 10.1088/00319155/49/18/012.PubMedGoogle Scholar
 Seo JK, Kwon O, Woo EJ: Magnetic resonance electrical impedance tomography (MREIT): conductivity and current density imaging. Journal of Physics: Conference Series. 2005, 12: 140155. 10.1088/17426596/12/1/014.Google Scholar
 Ozdemir M, Eyuboglu M, Ozbek O: Equipotential projectionbased magnetic resonance electrical impedance tomography and experimental realization. Physics in Medicine and Biology. 2004, 49: 47654783. 10.1088/00319155/49/20/008.PubMedGoogle Scholar
 Babiloni F, Mattia D, Babiloni C, Astolfi L, Salinari S, Basilisco A, Rossini PM, Marciani MG, Cincotti F: Multimodal integration of EEG, MEG and fMRI data for the solution of the neuroimage puzzle. Magn Reson Imaging. 2004, 22 (10): 14711476. 10.1016/j.mri.2004.10.007.PubMedGoogle Scholar
 Leahy RM, Mosher J, Spencer MC, Huang MX, Lewine J: A study of dipole localization accuracy for MEG and EEG using a human skull phantom. Electroencephalography and clinical Neurophysiology. 1998, 107 (2): 159173. 10.1016/S00134694(98)000571.PubMedGoogle Scholar
 Chang N, Gulrajani R, Gotman J: Dipole localization using simulated intracerebral EEG. Clin Neurophysiol. 2005, 116 (11): 27072716. 10.1016/j.clinph.2005.07.002.PubMedGoogle Scholar
 De Clercq W, Vergult A, Vanrumste B, Paesschen WV, Huffel SV: Canonical correlation analysis applied to remove muscle artifacts from the electroenceaphalogram. IEEE Transactions on Biomedical Engineering. 2006, 53 (12): 25832587. 10.1109/TBME.2006.879459.PubMedGoogle Scholar
 Zhukov L, Weinstein D, Johnson C: Independent component analysis for EEG source localization. Engineering in Medicine and Biology Magazine. 2000, IEEE, 19 (3): 8796. 10.1109/51.844386.Google Scholar
 Makeig S, Debener S, Onton J, Delorme A: Mining eventrelated brain dynamics. Trends in Cognitive Science. 2004, 8 (5): 204210. 10.1016/j.tics.2004.03.008.Google Scholar
 Astolfi L, Cincotti F, Mattia D, de Vico Fallani F, Salinari S, Ursino M, Zavaglia M, Marciani M, Babiloni F: Estimation of the cortical connectivity patterns during the intention of limb movements. Engineering in Medicine and Biology Magazine. 2006, IEEE, 25 (4): 3238. 10.1109/MEMB.2006.1657785.Google Scholar
 Gotman J, Kobayashi E, Bagshaw AP, Bnar CG, Dubeau F: Combining EEG and fMRI: a multimodal tool for epilepsy research. J Magn Reson Imaging. 2006, 23 (6): 906920. 10.1002/jmri.20577.PubMedGoogle Scholar
 Bénar CG, Grova C, Kobayashi E, Bagshaw AP, Aghakhani Y, Dubeau F, Gotman J: EEGfMRI of epileptic spikes: concordance with EEG source localization and intracranial EEG. Neuroimage. 2006, 30 (4): 11611170. 10.1016/j.neuroimage.2005.11.008.PubMedGoogle Scholar
 Pruis GW, Gilding BH, Peters MJ: A comparison of different numerical methods for solving the forward problem in EEG and MEG. Physiol Meas. 1993, 14 (Suppl 4A): A1A9. 10.1088/09673334/14/4A/001.PubMedGoogle Scholar
 Babuska I, Guo B: The hp version of the finite element method for problems with nonhomogeneous essential boundary condition. Computer Methods in Applied Mechanics and Engineering. 1989, 74: 128. 10.1016/00457825(89)900832.Google Scholar
 Melenk J, Schwab C: hpFEM for reactiondiffusion equations I: robust exponential convergence. SIAM journal on numerical analysis. 1998, 35: 15201557. 10.1137/S0036142997317602.Google Scholar
 Vejchodsky T, Solin P, Zitka M: Modular hpFEM System HERMES and its application to the Maxwell Equations. Mathematical Computations and Simulations. 2005,Google Scholar
 Wolters C, Köstler H, Möller C, Härdtlein J, Anwander A: Numerical approaches for dipole modeling in finite element method based source analysis. International Congres Series. 2007Google Scholar
Copyright
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.