Review on solving the forward problem in EEG source analysis

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 post-synaptic 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 three-shell 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 over-relaxation, 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.

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 so-called 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 event-related 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 state-of-the-art 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 three-shell 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 transfer-coefficients 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 over-relaxation (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 post-synaptic potentials and inhibitory post-synaptic 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¹⁰ 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 Goldman-Hodgkin-Katz 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 neurotrans-mitter. 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].

Figure 1 illustrates the excitatory and inhibitory postsynaptic potentials. It also shows the generation of an action potential. Further readings on the electrophysiology of neurons can be found in [2, 6].

Figure 1

Excitatory and inhibitory post synaptic potentials. An illustration of the action potentials and post synaptic potentials measured at different locations at the neuron. On the left a neuron is displayed and three probes are drawn at the location where the potential is measured. The above picture on the right shows the incoming exitatory action potentials measured at the probe at the top, at the probe in the middle the incoming inhibitory action potential is measured and shown. The neuron processes the incoming potentials: the excitatory action potentials are transformed into excitatory post synaptic potentials, the inhibitory action potentials are transformed into inhibitory post synaptic potentials. When two excitatory post synaptic potentials occur in a small time frame, the neuron fires. This is shown at the bottom figure. The dotted line shows the EPSP, in case there was no second excitatory action potential following. From [2].

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.

The following is focused on excitatory synapses and EPSP, located at the apical dendrites of a pyramidal cell. The neurotransmitter in the excitatory synapses causes an influx of positive ions at the postsynaptic membrane as illustrated in figure 2(a) and depolarizes the local cell membrane. This causes a lack of extracellular positive ions at the apical dendrites of the postsynaptic neuron. A redistribution of positively charged ions also takes place at the intracellular side. These ions flow from the apical dendrite to the cell body and depolarize the membrane potentials at the cell body. Subsequently positive charged ions become available at the extracellular side at the cell body and basal dendrites.

Figure 2

equivalent circuit for a neuron. An excitatory post synaptic potential, an simplified equivalent circuit for a neuron, and a resistive network for the extracellular environment. A neuron with an excitatory synapse at the apical dendrite is presented in (a). From [2]. A simplified equivalent circuit is depicted in (b). The extracellular environment can be represented with a resistive network as illustrated in (c).

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.

Quasi-static 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 quasi-static 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 current density is a vector field and can be represented by J(x, y, z). The unit of the current density is A/m². The divergence of a vector field J is defined as follows:

$$\nabla \cdot J=\underset{G\to 0}{\lim }\frac{1}{G}{\displaystyle {\oint }_{\partial G}JdS}$$

(1)

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³ 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₁(x₁, y₁, z₁) 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₁). 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.

For the third case a small volume around the current source at position r₂(x₂, y₂, z₂) is constructed. The current source represents the injection of positively charged ions at the cell body of the pyramidal cell. The current source density equals Iδ(r - r₂). Figure 3 represents the current density vectors for a current source and current sink configuration. Furthermore, three boxes are presented corresponding with the three cases discussed above.

Figure 3

The current density and equipotential lines in the vicinity of a dipole. The current density and equipotential lines in the vicinity of a current source and current sink is depicted. Equipotential lines are also given. Boxes are illustrated which represent the volumes G.

Uniting the three cases given above, one obtains:

∇·J = Iδ(r - r₂) - Iδ(r - r₁).

Ohm's law, the potential field and anisotropic/isotropic conductivities

The relationship between the current density J in A/m² and the electric field E in V/m is given by Ohm's law:

J = σ E

with σ(r) ∈ ℝ^3×3 being the position dependent conductivity tensor given by:

$$\sigma =\left[\begin{array}{ccc}{\sigma }_{11}& {\sigma }_{12}& {\sigma }_{13}\\ {\sigma }_{12}& {\sigma }_{22}& {\sigma }_{23}\\ {\sigma }_{13}& {\sigma }_{23}& {\sigma }_{33}\end{array}\right],$$

(5)

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).

At the skull, for example, the conductivity tangential to the surface is 10 times [16] the conductivity perpendicular to the surface (see figure 4(a)). The rationale for this is that the skull consists of 3 layers: a spongiform layer between two hard layers. Water, and also ionized particles, can move easily through the spongiform layer, but not through the hard layers [17]. Wolters et al. state that skull anisotropy has a smearing effect on the forward potential computation. The deeper a source lies, the more it is surrounded by anisotropic tissue, the larger the influence of the anisotropy on the resulting electric field. Therefore, the presence of anisotropic conducting tissues compromises the forward potential computation and as a consequence, the inverse problem [18].

Figure 4

Anisotropic conductivity of the brain tissues. The anisotropic properties of the conductivity of skull and white matter tissues. The anisotropic properties of the conductivity of skull and white matter tissues. (a) The skull consists of 3 layers: a spongiform layer between two hard layers. The conductivity tangentially to the skull surface is 10 times larger than the radial conductivity. (b) White matter consist of axons, grouped in bundles. The conductivity along the nerve bundle is 9 times larger than perpendicular to the nerve bundle.

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 (DT-MRI) [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.

In the grey matter, scalp and cerebro-spinal fluid (CSF) the conductivity is equal in all directions. Thus the place dependent conductivity tensor becomes a place dependent scalar σ, a so-called isotropic conducting tissue. The conductivity of CSF is quite accurately known to be 1.79 S/m [26]. In the following we will focus on the conductivity of the skull and soft tissues. Some typical values of conductivities can be found in table 1.

Table 1 The reference values of the absolute and relative conductivity of the compartments incorporated in the human head.

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 skull-to-soft 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, skull-to-soft 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 skull-to-soft tissue conductivity ratio was estimated to be around 20–25 [30, 31]. However in [30], it was shown that the skull-to-soft 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 skull-to-soft 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 skull-to-soft 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 soft-tissue 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 quasi-static 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₂) + Iδ(r - r₁).

In the Cartesian coordinate system equation (8) becomes for isotropic conductivities:

$$\begin{array}{ll}\frac{\partial }{\partial x}(\sigma \frac{\partial V}{\partial x})+\frac{\partial }{\partial y}(\sigma \frac{\partial V}{\partial y})+\frac{\partial }{\partial z}(\sigma \frac{\partial V}{\partial z})=\hfill & -I\delta (x-x_2)\delta (y-y_2)\delta (z-z_2)\hfill \\ +I\delta (x-x_1)\delta (y-y_1)\delta (z-z_1)\hfill \end{array}$$

(9)

and for anisotropic conductivities:

$$\begin{array}{c}{\sigma }_{11}\frac{{\partial }^{2}V}{\partial x^2}+{\sigma }_{22}\frac{{\partial }^{2}V}{\partial y^2}+{\sigma }_{33}\frac{{\partial }^{2}V}{\partial z^2}+2\left({\sigma }_{12}\frac{{\partial }^{2}V}{\partial x\partial y}+{\sigma }_{13}\frac{{\partial }^{2}V}{\partial x\partial z}+{\sigma }_{23}\frac{{\partial }^{2}V}{\partial y\partial z}\right)\\ +\left(\frac{\partial {\sigma }_{11}}{\partial x}+\frac{\partial {\sigma }_{12}}{\partial y}+\frac{\partial {\sigma }_{13}}{\partial z}\right)\frac{\partial V}{\partial x}+\left(\frac{\partial {\sigma }_{12}}{\partial x}+\frac{\partial {\sigma }_{22}}{\partial y}+\frac{\partial {\sigma }_{23}}{\partial z}\right)\frac{\partial V}{\partial y}+\left(\frac{\partial {\sigma }_{13}}{\partial x}+\frac{\partial {\sigma }_{23}}{\partial y}+\frac{\partial {\sigma }_{33}}{\partial z}\right)\frac{\partial V}{\partial z}=\\ -I\delta (x-x_2)\delta (y-y_2)\delta (z-z_2)+I\delta (x-x_1)\delta (y-y_1)\delta (z-z_1).\end{array}$$

(10)

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

At the interface between two compartments, two boundary conditions are found. Figure 5 illustrates such an interface. A first condition is based on the inability to pile up charge at the interface. All charge leaving one compartment through the interface must enter the other compartment. In other words, all current (charge per second) leaving a compartment with conductivity σ₁ through the interface enters the neighboring compartment with conductivity σ₂:

Figure 5

The boundary between two compartments. The boundary between two compartments. The boundary between two compartments, with conductivity σ₁ and σ₂. The normal vector e _n to the interface is also shown.

$$\begin{array}{c}{J}_1\cdot e_n=J_2\cdot e_n,\\ ({\sigma }_1\nabla V_1)\cdot e_n=({\sigma }_2\nabla V_2)\cdot e_n,\end{array}$$

(11)

where e _n is the normal component on the interface.

In particular no current can be injected into the air outside the human head due to the very low conductivity of the air. Therefore the current density at the surface of the head reads:

$$\begin{array}{c}{J}_1\cdot e_n=0,\\ ({\sigma }_1\cdot \nabla V_1)\cdot e_n=0.\end{array}$$

(12)

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₁ = V₂.

This equation represents the Dirichlet boundary condition.

The current dipole

Current source and current sink inject and remove the same amount of current I and they represent an active pyramidal cell at microscopic level. They can be modeled as a current dipole as illustrated in figure 6(a). The position parameter r _dip of the dipole is typically chosen half way between the two monopoles.

Figure 6

The dipole parameters. (a) The dipole parameters for a given current source and current sink configuration. (b) The dipole as a vector consisting of 6 parameters. 3 parameters are needed for the location of the dipole. 3 other parameters are needed for the vector components of the dipole. These vector components can also be transformed into spherical components: an azimuth, elevation and magnitude of the 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).

For N electrodes and p dipoles:

$$V=\left[\begin{array}{c}V(r_1)\\ ⋮\\ V(r_N)\end{array}\right]=\left[\begin{array}{ccc}g(r_1,r_{di{p}_1},e_{d_1})& \cdots & g(r_1,r_{di{p}_p},e_{d_p})\\ ⋮& \ddots & ⋮\\ g(r_N,r_{di{p}_1},e_{d_1})& \cdots & g(r_N,r_{di{p}_p},e_{d_p})\end{array}\right]\left[\begin{array}{c}{d}_1\\ ⋮\\ d_p\end{array}\right]=G(\{r_j,r_{di{p}_i},e_{d_i}\})\left[\begin{array}{c}{d}_1\\ ⋮\\ d_p\end{array}\right]$$

where i = 1,...,p and j = 1,...,N. Here V is a column vector.

For N electrodes, p dipoles and T discrete time samples:

$$\begin{array}{c}V=\left[\begin{array}{ccc}V(r_1,1)& \cdots & V(r_1,T)\\ ⋮& \ddots & ⋮\\ V(r_N,1)& \cdots & V(r_N,T)\end{array}\right]\\ =G(\{r_j,r_{di{p}_i},e_{d_i}\})\left[\begin{array}{ccc}{d}_{1,1}& \cdots & d_{1,T}\\ ⋮& \ddots & ⋮\\ d_{p,1}& \cdots & d_{p,T}\end{array}\right]=G(\{r_j,r_{di{p}_i},e_{d_i}\})D\end{array}$$

(17)

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 zeroth-order and first-order 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

The potential field generated by a current dipole with dipole moment d = d e _d at a position r _dip in an infinite conductor with conductivity σ, is introduced. The potential field is given by:

$$V(r,r_{dip},d)=\frac{d\cdot (r-r_{dip})}{4\pi \sigma \left|\right|r-r_{dip}|{|}^3},$$

(19)

with r being the position where the potential is calculated. Assume that the dipole is located in the origin of the Cartesian coordinate system and oriented along the z-axis. Then it can be written:

$$V(r,0,d{e}_z)=\frac{d\cos \theta }{4\pi \sigma r^2},$$

(20)

where θ represents the angle between the z-axis and r and r = ||r||. An illustration of the electrical potential field caused by dipole is shown in figure 7.

Figure 7

The equipotential lines of a dipole. The equipotential lines of a dipole oriented along the z-axis. The numbers correspond to the level of intensity of the potential field generated of the dipole. The zero line divides the dipole field into two parts: a positive one and a negative one.

Equation (20) shows that a dipole field attenuates with 1/r². 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

The first volume conductor models of the human head consisted of a homogeneous sphere [43]. However it was soon noticed that the skull tissue had a conductivity which was significantly lower than the conductivity of scalp and brain tissue. Therefore the volume conductor model of the head needed further refinement and a three-shell concentric spherical head model was introduced. In this model, the inner sphere represents the brain, the intermediate layer represents the skull and the outer layer represents the scalp. For this geometry a semi-analytical solution of Poisson's equation exists. The derivation is based on [44, 45]. Consider a dipole located on the z-axis and a scalp point P, located in the xz-plane, as illustrated in figure 8. The dipole components located in the xz-plane i.e. d _r the radial component and d _t the tangential component, are also shown in figure 8. The component orthogonal to the xz-plane, does not contribute to the potential at scalp point P due to the fact that the zero potential plane of this component traverses P. The potential V at scalp point P for the proposed dipole is given by:

$$V=\frac{1}{4\pi S{R}^2}{\displaystyle \sum _{i=1}^{\infty }\frac{X{(2i+1)}^3}{g_i(i+1)i}{b}^{i-1}[i{d}_r{P}_i(\cos \theta )+d_t{P}_i^1(\cos \theta )],}$$

(21)

Figure 8

The three-shell concentric spherical head model. The dipole is located on the z-axis and the potential is measured at scalp point P located in the xz-plane.

with g _i given by:

$$g_i=[(i+1)X+i][\frac{iX}{i+1}+1]+(1-X)[(i+1)X+i](f_1^{i_1}-f_2^{i_1})-i{(1-X)}^2{(f_1/f_2)}^{i_1}.$$

(22)

Where:

d _r is the radial component (3 × 1-vector in meters),

d _t is the tangential component (3 × 1-vector 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₁ equals 2i + 1,

r₁ is the radius of the inner shell (in meters),

r₂ is the radius of the middle shell (in meters),

f₁ equals r₁/R (unitless) and

f₂ equals r₂/R (unitless).

Equation (21) gives the scalp potentials generated by a dipole located on the z-axis, 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 semi-analytical 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 three-shell spherical head model, such as the Berg approximation [53], in which a single-sphere model is used to approximate a three- (or four-) layer sphere, have also been used to improve further the computational efficiency of multi-layer 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. So-called real head models are becoming much more common in the literature, in conjunction with either boundary-element, finite-element, or finite-difference methods. However, the computational requirements for a realistic head model are higher than that for a multi-layer sphere.

An approach which is situated between the spherical head model approaches and realistic ones is the sensor-fitted 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 non-conducting 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: brain-skull interface, skull-scalp 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.

The integral equations describing the potential V(r) at any point r in a piecewise volume conductor V were described in [61–63]:

$$V(r)=\frac{2{\sigma }_0}{{\sigma }_k^{-}+{\sigma }_k^{+}}{V}_0(r)+\frac{1}{2\pi }{\displaystyle \sum _{j=1}^R\frac{{\sigma }_j^{-}\_{\sigma }_j^{+}}{{\sigma }_k^{-}+{\sigma }_k^{+}}{\displaystyle {\int }_{r^{\prime }\in S_j}V(r^{\prime })\frac{r^{\prime }-r}{\left|\right|r^{\prime }-r|{|}^3}d{S}_j}},$$

(23)

where σ₀ corresponds to the medium in which the dipole source is located (the brain compartment) and V₀(r) is the potential at r for an infinite medium with conductivity σ₀ 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.

Each interface S _j is digitized in $N_{S_j}$ triangles, (see figure 9) and in each triangle centre the potentials are calculated using equation (23). The integral over the surface S _j is transformed into a summation of integrals over traingles on that surface. The potential values on surface S _j can be written as

$$V(r)=\frac{2{\sigma }_0}{{\sigma }_r^{-}+{\sigma }_r^{+}}{V}_0(r)+\frac{1}{2\pi }{\displaystyle \sum _{k=1}^R\frac{{\sigma }_k^{-}\_{\sigma }_k^{+}}{{\sigma }_r^{-}+{\sigma }_r^{+}}{\displaystyle \sum _{j=1}^{N_{S_k}}{\displaystyle {\int }_{{\Delta }_{S_k,j}}V(r^{\prime })\frac{r^{\prime }-r}{\left|\right|r^{\prime }-r|{|}^3}d{S}_k}}},$$

(24)

Figure 9

Example mesh of the human head used in BEM. Triangulated surfaces of the brain, skull and scalp compartment used in BEM. The surfaces indicate the different interfaces of the human head: air-scalp, scalp-skull and skull-brain.

where the integral is over ${\Delta }_{S_j,k}$, the j-th triangle on the surface S _j , R is the number of interfaces in the volume. An exact solution of the integral is generally not possible, therefore an approximated solution ${\tilde{V}}^k(r)$ on surface S _k may be defined as a linear combination of $N_{S_k}$ simple basis functions

$$\widetilde{V^k}(r)={\displaystyle \sum _{i=1}^{N_{S_k}}{V}_i^k{h}_i(r).}$$

(25)

The coefficients $V_i^k$ represent unknowns on surface S _k whose values are determined by constraining $\tilde{V}(x)$ to satisfy (24) at discrete points, also known as collocation points. Moreover, equation (24) can be rewritten as

$$V(r)=\frac{2{\sigma }_0}{{\sigma }_r^{-}+{\sigma }_r^{+}}{V}_0(r)+\frac{1}{2\pi }{\displaystyle \sum _{k=1}^R\frac{{\sigma }_k^{-}\_{\sigma }_k^{+}}{{\sigma }_r^{-}+{\sigma }_r^{+}}{\displaystyle \sum _{j=1}^{N_{S_k}}{\displaystyle \sum _{i=1}^{N_{S_k}}{V}_i^k}{\displaystyle {\int }_{{\Delta }_{S_k,j}}{h}_i(r)\frac{r^{\prime }-r}{\left|\right|r^{\prime }-r|{|}^3}d{S}_k}}}.$$

(26)

This equation can be transformed into a set of linear equations:

V = BV + V₀,

where V and V₀ 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.

Determination of the elements of the matrix B is computationally intensive and there exist different approaches for their computation. The integral in equation (23) is also often called the solid angle [62, 64, 65]. The basis functions h _i (r) can be defined in several ways. The "constant-potential" approach for triangular elements uses basis functions defined by

$$h_i(r)=\{\begin{array}{cc}1& r\in {\Delta }_i\\ 0& r\notin {\Delta }_i\end{array}$$

(28)

where Δ _i denotes the ith planar triangle on the tesselated surface. The collocation points are typically the centroids of the surface elements and the unknown potentials V are the potentials at each triangle [66]. The "linear potential" approach uses basis functions defined by

$$h_i(r)=\{\begin{array}{ll}\frac{[r{r}_j{r}_k]}{[r_i{r}_j{r}_k]}\hfill & r\in {\Delta }_{i(jk)}\hfill \\ 0\hfill & r\notin {\Delta }_{i(jk)}\hfill \end{array}$$

(29)

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 higher-order 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].

Barnard et al. [64] showed that the potentials in equation (27) are only defined up to an additive constant. Hence, equation (27) has no unique solution. This ambiguity can be removed by deflation, which means that B must be replaced by

$$C=B-\frac{1}{N}e{e}^T,$$

(30)

where e is a vector with all its N (the total number of unknowns) components equal to one. The deflated equation

V = CV + V₀,

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^-1V₀.

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'' is the potential on surface S _brain when the head is a homogeneous brain region, thus omitting the skull and scalp compartments. V' is the correction term. When V is written like above, equation (32) can be written as

$$\begin{array}{c}{V}^{\prime }+{V^{\prime }}^{\prime }=A^{-1}{V}_0\\ V^{\prime }=A^{-1}(V_0-A{V^{\prime }}^{\prime })\\ V^{\prime }=A^{-1}{V^{\prime }}_0.\end{array}$$

(34)

Because V'' is zero on the interfaces S _scalp and S _skull , V' contains the potential values on the outer surface. The IPA is based on the more accurate solution of the right-hand side term ${V^{\prime }}_0$. An accurate solution can be obtained by setting ${V^{\prime }}_0$ to the following

$${V^{\prime }}_0=\left(\begin{array}{c}{V^{\prime }}_0^1\\ {V^{\prime }}_0^2\\ {V^{\prime }}_0^3\end{array}\right)=\left(\begin{array}{c}\beta V_0^1\\ \beta V_0^2\\ \beta V_0^3-\frac{2\beta }{\beta +1}{{V^{\prime }}^{\prime }}_0^3\end{array}\right),$$

(35)

where $V_0^1$, $V_0^2$ and $V_0^3$ are the potentials at respectively the brain-skull surface, the skull-scalp 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 multi-sphere models by Gençer and Akahn-Acar [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² 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² 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. Non-nested compartments are compartments that are not part of the brain, but part of the head (such as eyes, sinuses,...) [76].

The finite element method

Another method to solve Poisson's equation in a realistic head model is the finite element method (FEM). The Galerkin approach [77] is used to equation (7) with boundary conditions (11), (12), (13). First, equation (7) is multiplied with a test function φ and then integrated over the volume G representing the entire head. Using Green's first identity for integration:

$${\displaystyle {\int }_G\nabla \phi \cdot (\sigma \nabla V)dG=}{\displaystyle {\int }_{\partial G}\phi \sigma \nabla V\cdot dS}-{\displaystyle {\int }_G\phi (\nabla \sigma \nabla V)dG},$$

(36)

in combination with the boundary conditions (12), yields the 'weak formulation' of the forward problem:

$$-{\displaystyle {\int }_G\nabla \phi \cdot (\sigma \nabla V)dG}={\displaystyle {\int }_G\phi I_{m}dG}.$$

(37)

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 , φ)

The entire 3D volume conductor is digitized in small elements. Figure 10 illustrates a 2D volume conductor digitized with triangles.

Figure 10

Example mesh in 2D used in FEM. A digitization of the 2D coronal slice of the head. The 2D elements are the triangles.

The computational points ${\left\{V_i\right\}}_{i=1}^n$ can be identified with the vertices of the elements (n is the number of vertices). The unknown potential V(x, y, z) is given by

$$V(x,y,z)={\displaystyle \sum _{i=1}^n{V}_i{\phi }_i(x,y,z)},$$

(39)

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 non-zero is limited to adjacent elements. Moreover, the basis functions span a space of piecewise polynomial functions.

Furthermore, they have the property that they are each equal to unity at the corresponding computational point and equal to zero at all other computational points. Substituting (39) in equation (38) produces n equations in n unknowns V = [V₁...V _n ]^T∈ ℝ^n×1:

$$a({\displaystyle \sum _{i=1}^n{V}_i{\phi }_i,{\phi }_j})=(I_m,{\phi }_j),$$

(40)

$${\displaystyle \sum _{i=1}^{n}a({\phi }_i,{\phi }_j})V_i=(I_m,{\phi }_j)$$

(41)

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 finite-element 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 FEM-mesh. 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 4-layer 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 equivalent-current 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 finite-element 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 equivalent-dipole inverse solution (EDS) in that it provides source-sink 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)

The differential equation (9) with boundary conditions (11), (12), (13) is transformed into a linear equation utilizing the 'box integration' scheme [86] for the cell-centered iFDM. Consider a typical node P in a cubic grid with internode spacing h. The six neighbouring nodes are Q _i (i = 1,...,6) as illustrated in figure 11.

Figure 11

The computation stencil used in FDM. A typical node P in an FDM grid with its neighbours Q _i (i = 1⋯6). The volume G₀ is given by the box.

Introducing α _i and α₀ as,

$$\begin{array}{lll}{\alpha }_i\hfill & =\hfill & 2h\frac{{\sigma }_0{\sigma }_i}{{\sigma }_0+{\sigma }_i}\hfill \\ {\alpha }_0\hfill & =\hfill & {\displaystyle \sum _{i=1}^6{\alpha }_i},\hfill \end{array}$$

(43)

a finite difference approximation of (9) is obtained:

$${\displaystyle \sum _{i=1}^6{\alpha }_i{V}_{Q_i}}-{\alpha }_0{V}_P=I_P,$$

(44)

with

$$I_P={\displaystyle \int {\displaystyle \int {\displaystyle {\int }_{G_0}-I\delta (x-x_2)\delta (y-y_2)\delta (z-z_2)+I\delta (x-x_1)\delta (y-y_1)\delta (z-z_1)dxdydz}}.}$$

(45)

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×nare 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×nhas at most six off-diagonal 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.

If the local coordinate system coincides with the axes representing the principal directions of the anisotropy, then the conductivity tensor at an element can be written as a diagonal matrix. The diagonal elements represent the conductivities in the orthogonal directions. The matrix representation has to be transformed to the global Cartesian system of the head, the same for all elements. A rotation matrix is then required to transform the principal directions to a conductivity tensor in the Cartesian coordinate system. In the local coordinate system the conductivity tensor at an element j can be written as follows:

$${\sigma }^j=\left(\begin{array}{lll}{\sigma }_1^i\hfill & 0\hfill & 0\hfill \\ 0\hfill & {\sigma }_2^i\hfill & 0\hfill \\ 0\hfill & 0\hfill & {\sigma }_3^i\hfill \end{array}\right),$$

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 ^Tdenotes 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].

A typical node 0 in the grid represents the intersection of eight neighbouring cubic elements, as shown in figure 12. The finite difference formulation of equation 10 at node 0, derived from Saleheen's method. From the 26 neigbouring nodes at node 0, the formulation uses the 18 nearest neighbours, with rectilinear distance ≤ 2:

$${\displaystyle \sum _{i=1}^{18}{a}_i{V}_i}-\left({\displaystyle \sum _{i=1}^{18}{a}_i}\right)V_0=I,$$

(47)

where V _i is the discrete potential value at node i. a _i are the coefficients depending on the conductivity tensor of the elements and the internode distance. These coefficients are given in [96].

Figure 12

The computation stencil used in FDM if anistropic conductivities are incorporated. The potential at node 0 can be written as a linear combination of 18 neighbouring nodes in the FDM scheme. For each node we obtain an equation, which can be put into a linear system Ax = b.

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

The three methods BEM, FEM and FDM can all be used to solve the forward problem of EEG source analysis in a realistic head model. A summary of the comparison between the BEM, FEM and FDM is given below and in table 2.

Table 2 A comparison of the different methods for solving Poisson's equation in a realistic head model is presented.

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 ₀ 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 right-hand side term, utilizing iterative solvers such as the successive over-relaxation 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^-1Ax = M^-1b, 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, Gauss-Seidel, Successive Over-Relaxation (SOR) and Symmetric Successive Over-Relaxation (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 transfer-coefficients 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 transfer-coefficients 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

Consider a resistor circuit, with two clamps AB and r _x as illustrated in figure 13. The clamp AB represents a pair of scalp electrodes. The clamp r _x is located in the brain region.

Figure 13

Reciprocity. A resistor network where a current source is introduced in the brain and the a potential difference is measured at an electrode pair, and visa versa: (a) a current source $I_{r_x}$ is introduced and the potential U _AB is measured, and (b) a current source I _AB is introduced and a potential $V_{r_x}$ is measured.

First a current $I_{r_x}$ at clamp r _x is introduced. This source will generate a potential U _AB ($I_{r_x}$) at AB as illustrated in figure 13(a). Next, current I _AB at AB is set. This will give rise to a potential difference $V_{r_x}$(I _AB ) at r _x illustrated in figure 13(b). The reciprocity theorem in circuit analysis states:

$$U_{AB}{I}_{AB}=V_{r_x}{I}_{r_x}.$$

(48)

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×nand has the following properties: it is sparse, symmetric and regular. One can write:

A·V = I,

with V = [V₁...V _n ]^T∈ ℝ^n×1 and I = [I₁...I _n ]^T∈ ℝ^n×1 and with ^Tthe transpose operator. The desired potential difference V _k - V _l between nodes k and l can be obtained for a current source I _f at node f and a current sink I _g at node g with I _f = -I _g . All other sources are zero. Cramer's solution for a linear system then becomes:

$$V_k=\frac{I_f[{(-1)}^{k+f+1}{A}_{fk}-{(-1)}^{k+g+1}{A}_{gk}]}{\det A},$$

(49)

$$V_l=\frac{I_f[{(-1)}^{l+f+1}{A}_{fl}-{(-1)}^{l+g+1}{A}_{gl}]}{\det A},$$

(50)

with A_*∘ the minor for row * and column ∘.

On the other hand the potential V _f and V _g for a current source I _k and current sink I _l with I _k = -I _l , are:

$$V_f=\frac{I_k[{(-1)}^{f+k+1}{A}_{kf}-{(-1)}^{f+l+1}{A}_{lf}]}{\det A},$$

(51)

$$V_g=\frac{I_k[{(-1)}^{g+k+1}{A}_{kg}-{(-1)}^{g+l+1}{A}_{lg}]}{\det A}.$$

(52)

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

Considering equation (48), a dipole can be represented as two current monopoles, a current source and sink, providing $I_{r_x}$ and -$I_{r_x}$, separated by a distance 2h. The dipole is oriented from the negative to the positive current monopole and is assumed to be along the x-axis of the resistor network with node spacing h. The magnitude of the dipole moment is then 2h$I_{r_x}$. The centre r of the two monopoles can then be seen as the dipole position. The scalp electrodes are located sufficiently far from the sources compared with the distance 2h between the sources so that we can assume a dipole field. Equation (48) can be rewritten as:

$$U_{AB}=\frac{V_{r_x}{I}_{r_x}}{I_{AB}}.$$

(53)

The forward problem in EEG source analysis gives the potential U _AB for a current dipole located at r and oriented along the x-axis. Rewriting equation(53) with d _x = 2h$I_{r_x}$ and

$$\frac{\partial V}{\partial x}\approx \frac{[V_{I_{AB}}(r+h{e}_x)-V_{I_{AB}}(r-h{e}_x)]}{2h},$$

gives:

$$U_{AB}=\frac{d_x\frac{\partial V}{\partial x}}{I_{AB}}.$$

(54)

In a similar way, U _AB can be calculated for a dipole located at r oriented along the y-axis and the z-axis.

Consider a dipole at position r and with dipole components d = (d _x , d _y , d _z )^T∈ ℝ^3×1. The potential U _AB reads:

$$U_{AB}(r,d)=\frac{d^T\cdot \nabla V(r)}{I_{AB}},$$

(55)

with ∇V(r) = (∂V(r)/∂x, ∂V(r)/∂y, ∂V(r)/∂z)^T∈ ℝ^3×1.

The flowchart in figure 14 shows the consecutive steps that are necessary in the reciprocity approach in conjunction with FDM.

Figure 14

The consecutive steps when applying reciprocity in conjunction with FDM. A scheme showing the consecutive steps that have to be followed when applying reciprocity in conjunction with FDM. First a current dipole I _AB is set on the electrode pair AB. Using FDM the potential field is calculated in each point V(ih, jh, kh). With the dipole parameters and the potential field, the reciprocity theorem can be applied. This results in a potential difference at the electrode pair AB.

• 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.

Figure 15

Lead field between two electrodes. The current density J = σ∇V and the equipotential lines are illustrated when introducing a current I _AB at electrode A and removing the same amount at electrode B.

• 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 tri-linear 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].

For FDM in isotropic or anisotropic media, it can be shown from equation (44) that the sum of all entries in a row/column of A equals zero (see equations 44 and 47). Therefore, the vector e = [1,...,1]^Tis an eigenvector with associated eigenvalue 0. The matrix (A) of the FDM in both isotropic and anisotropic media has rank n - 1, with n the number of unknowns, and the eigenspace of the eigenvalue 0 is of dimension one. Note that for a singular problem to have a solution at all, the right-hand side b must be consistent, i.e. b ∈ Range(A), Range(A) being the range of A and defined as the number of independent vectors in A. The kernel of A, Kernel(A), is the set of vectors a that if multiplied by A returns zero. In this case the problem Ax = b possesses an infinite set of solutions. An iterative method that converges from each initial guess towards an element of this solution set is said to be semi-convergent [100]. In our case, A is symmetric, thus Range(A) = Kernel(A)^⊥ where ^⊥ stands for the orthogonal complement. Since Kernel(A) is spanned by the vector e = (1,...,1)^Tcontaining only ones as entries, a vector v lies in Range(A) if and only if

$$e^{T}v={\displaystyle \sum _{k=1}^n{v}_k=0}.$$

(56)

From equations (44) and (47) it is easy to see that the right-hand 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 right-hand 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 non-zero off-diagonal elements. This leads to a very small ratio of non-zero to overall entries resulting in a very sparse matrix.

Iterative solvers

The following methods will be discussed:

• Successive over-relaxation (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 semi-definite case. In the case of a consistent right-hand side, semi-convergence 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 over-relaxation

The SOR method is a representative of classical stationary methods. It is known to be a non-optimal 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_i1x₁ + ⋯ + a _ii x _i + ⋯ + a _in x _n = b _i ,   i = 1,...,n,

can be rewritten as

$$\begin{array}{cc}{x}_i=\frac{1}{a_{ii}}\left(b_i-{\displaystyle \sum _{j=1,j\ne i}^n{a}_{ij}{x}_j}\right).& i=1,\mathrm{...},n,\end{array}$$

Let x^(k) be an approximation to the solution after k iterations. The SOR method updates the unknowns in the following fashion. To compute $x_i^{(k+1)}$ first an intermediate value

$${\overline{x}}_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-{\displaystyle \sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}}-{\displaystyle \sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}}\right)$$

(57)

is determined. Here new values of x^(k+1) are used as soon as they are available. The new approximation then becomes

$$x_i^{(k+1)}=\omega {\overline{x}}_i^{(k+1)}+(1-\omega )x_i^{(k)}=x_i^{(k)}+\omega \left({\overline{x}}_i^{(k+1)}-x_i^{(k)}\right).$$

(58)

The over-relaxation parameter ω is a weighting parameter used to put more weight onto the correction in order to improve convergence. According to the Young theorem, the optimal value for ω can be computed and can be shown to be equal to:

$${\omega }_{opt}=\frac{2}{1+\sqrt{1-\rho {(B)}^2}},$$

(59)

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].

The pseudocode for the SOR algorithm is given in figure 16.

Figure 16

The SOR method. Pseudo-code for the successive over-relaxation method. The instructions to be processed in a for-loop are indicated between the do and od.

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),

where d^(k) ∈ ℝⁿis a search direction and α^(k) is a scalar given by

$${\alpha }^{(k)}=\frac{{(r^{(k)})}^T{r}^{(k)}}{{(d^{(k)})}^{T}A{d}^{(k)}}.$$

The first search direction is just the residual of the initial guess d⁽⁰⁾ = r⁽⁰⁾, where the residual is defined by r^(k) = b - Ax^(k). The k-th search direction is computed from the previous one via

d^(k) = r^(k) + β^(k)d^(k-1),

with

$${\beta }^{(k)}=\frac{{(r^{(k)})}^T{r}^{(k)}}{{(r^{(k-1)})}^T{r}^{(k-1)}}.$$

The pseudocode of the CG method is given in figure 17 with M equal to the unit matrix. Let us turn our attention now to the use of symmetric successive over-relaxation (SSOR) as a preconditioner for the CG method.

Figure 17

The (P)CG method. Pseudo-code for the preconditioned conjugate gradient method. The instructions to be processed in a for-loop are indicated between the do and od.

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^-1A(E^-1)^T, $\tilde{x}$ = E^Tx and $\tilde{b}$ = E^-1b. 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^Tand with r the residual. Compare the pseudocode of the preconditioned CG in figure 17.

Algebraic multigrid

The last contestant is an algebraic multigrid method. Algebraic multigrid methods are known to be, in general, very efficient solvers for elliptic boundary value problems. The basic idea is the recursive application of a two-grid method. Here one splits the error into two components. These are typically referred to as rough and smooth, because in traditional applications they represent high- and low-frequency Fourier modes. The rough components are reduced in size on the original (fine) grid by applying a limited number of steps of some iteration scheme, such as Gauß-Seidel or SOR. This process is called smoothing, because in the remaining error the smooth components are dominant. For these a correction is computed on a coarser grid with a larger mesh size. The equation for this correction is then again solved by a two-grid approach, so that a hierarchy of grid levels is obtained. Figure 18 illustrates the pseudocode of one iteration of the algebraic multigrid method. Here $I_h^H$ represents the transfer function from a fine grid (h) to the next coarser grid (H) and $I_H^h$ from a coarse grid to the next finer grid. Furthermore c* represents the correction on grid '*' applied to update the solution on the next finer grid.

Figure 18

The Multigrid V-cycle. Pseudo-code of the Multigrid V-cylce. The instructions to be processed in a for-loop are indicated between the do and od.

The difficulty in algebraic multigrid is finding the proper components, that is coarsening strategies to derive a suitable grid hierarchy, operators for transferring functions between different grid levels, and smoothing iterations for the rough components. In the case of complex geometries and/or jumping coefficients this can be quite tedious. Therefore the idea of algebraic multigrid methods, as illustrated in [109, 110] is again attracting increased attention. Here a "grid hierarchy" and inter-grid transfer functions are derived automatically from the problem matrix. The pseudocode of the AMG method is given in figure 19.

Figure 19

The AMG method. Pseudo-code of the algebraic multigrid method. The instructions to be processed in a for-loop are indicated between the do and od.

Comparing the iterative solvers

For the iFDM, the following conclusions by [111] were drawn from comparing the iterative solvers: (a) the algebraic multigrid-based solver performed the task 1.8–3.5 times faster, platform depending, than the second-best contender, (b) there is no need to introduce a reference potential which forces a unique solution and (c) neither the grid- nor matrix-based implementation of the solvers consistently gives a smaller run time.

Wolters et al. [79] investigated the parallel implementation of iterative solvers. If the Jacobi-CG (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 AMG-CG (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 state-of-the-art methods applied to solve the forward problem. Multiple references to the work of authors active in this area were provided.

The post-synaptic 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 quasi-static 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 so-called forward problem. The first models used were three-shell 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 transfer-coefficients 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 right-hand 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 right-hand side. The successive over-relaxation method, the conjugate gradient method, the preconditioned conjugate gradient method, and the multigrid method were discussed. The last method was the most promising computation-time 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 bulk-conductivity 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 current-voltage 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/μ₀, μ₀ 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 high-frequent 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 long-term monitoring. Nevertheless, the combined use of EEG and MEG has shown to be beneficial for stimulus-locked 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 low-pass 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 signal-to-noise 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 extra-cerebral 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 h-version) FEM, while increasing the order of an element is known as the p-version of FEM. The generalisation of both variants is the hp-version of FEM (hp-FEM) 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 hp-FEM 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 right-hand 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.