2. The mathematical formulation of the EEG problem

Think of the following scenario. You are conducting a series of experimental or clinical studies, but out of curiosity and foremost for a better understanding of the way the system under consideration behaves, you desire to build a mathematical model interpreting as best as possible your measurements. But, where to begin? First of all, we will need a framework which is capable of, at least to a degree, explaining what happens and why. If such a framework does not exist, we have to formulate one. As mentioned in Section 1, the electrochemical activity of brain cells results in bioelectric sources which generate an electric field in the neighbourhood of the cells. This field varies generally in time. Consequently, electromagnetic phenomena materialize. Luckily for us, the framework interpreting this kind of phenomena and the starting point of our endeavour are Maxwell's equations, a set of four equations, namely

$$
\nabla \times \mathbf{E} + \frac{\partial \mathbf{B}}{\partial t} = \mathbf{0}, \ \nabla \times \mathbf{H} = \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t}, \ \nabla \cdot \mathbf{B} = \mathbf{0}, \ \nabla \cdot \mathbf{D} = \rho. \tag{1}
$$

The inverted delta present in Eq. (1) is called del, or nabla, and consist of a mathematical device named operator, a symbol indicating that an action must be performed on what follows. The algebraic operations of the dot (�) and cross (�) product between two vectors should not be confused with the elementary operation of multiplication.

Maxwell equations connect the electric fields E and B, the displacement field D and the magnetizing field H with their sources, that is, the charge density ρ and the current density J: These fields have direction and magnitude and must be represented as vector functions (bold capital letters). The above set of equations is also called Maxwell's macroscopic equations, and in order to apply them, a relation between D and E, as well as H and B must be specified. For materials without polarization and magnetization, they are

$$\mathbf{D} = \varepsilon \mathbf{E}, \ \mathbf{H} = \mu^{-1} \mathbf{B}, \tag{2}$$

where ε is the relative permittivity (dielectric constant) of the material and shows how strong the material influences the electric field E: Similarly, μ is the magnetic permeability of the material and provides the corresponding influence on the magnetic field B:

In what follows, we shall consider the following instance. For a finite medium, we introduce the characteristic dimension R, namely the smallest sphere with radius R which envelopes the medium under consideration. On the special occasion where the wavelength λ of the wave generated by the electromagnetic field is much larger than the characteristic dimension of the medium, namely λ>>R, then the corresponding time rates of change are very small. The latter observation leads to the quasi-static theory of Maxwell's equations, which take the form (see Ref. [18] for details)

$$
\nabla \times \mathbf{E} = \mathbf{0}, \ \nabla \times \mathbf{H} = \mathbf{J}, \ \nabla \cdot \mathbf{B} = 0, \ \nabla \cdot \mathbf{D} = \rho. \tag{3}
$$

It can be shown [19] that for a medium the size of the brain, R equals about 20 cm whereas λ about 400 m. Therefore, λ ffi 2000R and the application of Maxwell's quasi-static equations are justified. Replacing Eq. (2) into Eq. (3), we immediately find

$$
\nabla \times \mathbf{E} = \mathbf{0}, \ \nabla \times \mathbf{B} = \mu \mathbf{J}, \ \nabla \cdot \mathbf{B} = 0, \ \nabla \cdot \mathbf{D} = \varepsilon^{-1} \rho. \tag{4}
$$

Let us focus on the first of Eq. (4). The curl of a vector, that is, the operator ∇� captures the idea of how the vector field is circulating around a central axis. From this point of view, the electric field is irrotational. By a well-known identity of vector calculus, if the curl of a vector is zero, then the vector field in question can always be expressed as the gradient of a scalar field. In our case, the absence of circulation of E is caused by a continuously decreasing electric potential U along the direction of the electric field,

$$\mathbf{E} = -\nabla \mathbf{U}.\tag{5}$$

An inherit characteristic of linear systems is the principle of superposition. Here, the electric fields are superposable, meaning that the electric field generated by a number of charges can be expressed as the vector sum of the electric fields generated by each charge separately; it follows from Eq. (5) that the electric potentials are superposable as well. As a result, it is much easier to compute the provoked potential than the corresponding electric field.

At this point, the benefit of Eq. (5) is not clear yet. To show the usefulness of Eq. (5), we utilize a theorem of vector calculus stating that the divergence of the curl of a vector always vanishes. Applying the latter on the second of Eq. (3), we are left with

$$\nabla \cdot \mathbf{J} = 0 \tag{6}$$

and the current density of J is said to be solenoidal and expresses the steady-state condition that the charge density ρ is not changing in time. Moreover, an applied field in a resistive material, such as the brain, will induce a current of density J<sup>i</sup> , directly proportional to the applied field provided by a generalization of Ohm's law due to Kirchhoff (for a detailed historical analysis see Ref. [20]),

$$\mathbf{J}\_{\mathrm{i}} = \sigma \mathbf{E} \,\mathrm{\,} \tag{7}$$

where the proportionality constant σ is termed conductivity. Therefore, the total current is given as the sum of the primary current J<sup>p</sup> responsible for the electric field, and the induced, or secondary, current J<sup>i</sup> given by Eq. (8), as

$$\mathbf{J} = \mathbf{J}\_{\mathbf{p}} + \sigma \mathbf{E}.\tag{8}$$

Combining Eqs. (5), (6) and (8), we arrive at the following differential equation:

$$\nabla \cdot (\sigma \,\nabla U) = \nabla \cdot \mathbf{J}\_{\mathbb{P}'} \tag{9}$$

which must be satisfied by the electric potential U: Simple vector calculus shows that

$$\nabla \cdot (\sigma \,\nabla \mathcal{U}) = \nabla \mathcal{U} \cdot \nabla \sigma + \sigma \,\Delta \mathcal{U},\tag{10}$$

where the symbol Δ is called the Laplacian operator.

algebraic operations of the dot (�) and cross (�) product between two vectors should not be

Maxwell equations connect the electric fields E and B, the displacement field D and the magnetizing field H with their sources, that is, the charge density ρ and the current density J: These fields have direction and magnitude and must be represented as vector functions (bold capital letters). The above set of equations is also called Maxwell's macroscopic equations, and in order to apply them, a relation between D and E, as well as H and B must be specified. For

<sup>D</sup> <sup>¼</sup> <sup>ε</sup>E, <sup>H</sup> <sup>¼</sup> <sup>μ</sup>�<sup>1</sup>

material and provides the corresponding influence on the magnetic field B:

where ε is the relative permittivity (dielectric constant) of the material and shows how strong the material influences the electric field E: Similarly, μ is the magnetic permeability of the

In what follows, we shall consider the following instance. For a finite medium, we introduce the characteristic dimension R, namely the smallest sphere with radius R which envelopes the medium under consideration. On the special occasion where the wavelength λ of the wave generated by the electromagnetic field is much larger than the characteristic dimension of the medium, namely λ>>R, then the corresponding time rates of change are very small. The latter observation leads to the quasi-static theory of Maxwell's equations, which take the form (see Ref.

It can be shown [19] that for a medium the size of the brain, R equals about 20 cm whereas λ about 400 m. Therefore, λ ffi 2000R and the application of Maxwell's quasi-static equations are

<sup>∇</sup> � <sup>E</sup> <sup>¼</sup> <sup>0</sup>, <sup>∇</sup> � <sup>B</sup> <sup>¼</sup> <sup>μ</sup>J, <sup>∇</sup> � <sup>B</sup> <sup>¼</sup> <sup>0</sup>, <sup>∇</sup> � <sup>D</sup> <sup>¼</sup> <sup>ε</sup>�<sup>1</sup>

Let us focus on the first of Eq. (4). The curl of a vector, that is, the operator ∇� captures the idea of how the vector field is circulating around a central axis. From this point of view, the electric field is irrotational. By a well-known identity of vector calculus, if the curl of a vector is zero, then the vector field in question can always be expressed as the gradient of a scalar field. In our case, the absence of circulation of E is caused by a continuously decreasing electric potential U

An inherit characteristic of linear systems is the principle of superposition. Here, the electric fields are superposable, meaning that the electric field generated by a number of charges can be expressed as the vector sum of the electric fields generated by each charge separately; it

∇ � E ¼ 0, ∇ � H ¼ J, ∇ � B ¼ 0, ∇ � D ¼ ρ: ð3Þ

E ¼ �∇U: ð5Þ

B, ð2Þ

ρ: ð4Þ

confused with the elementary operation of multiplication.

materials without polarization and magnetization, they are

justified. Replacing Eq. (2) into Eq. (3), we immediately find

along the direction of the electric field,

[18] for details)

38 Electroencephalography

When the conductivity varies in space, that is, the medium under consideration is inhomogeneous, consisting of compartments which are not of the same material, the quantity ∇σ differs from zero. On the other hand, when the conductivity is constant in space (homogeneous, isotropic) or constant by direction (anisotropic), then the gradient of σ vanishes. In the latter case, the electric potential U is related to the primary current J<sup>p</sup> by Poisson's equation,

$$
\Delta \mathcal{U} = \sigma^{-1} \nabla \cdot \mathbf{J}\_{\mathbb{P}}.\tag{11}
$$

Note that the 'source term', namely the right-hand side of Eq. (11), is not the primary current per se, but the divergence of J<sup>p</sup> multiplied by 1 over σ. In view of Eq. (11), the forward EEG problem is now formulated as follows. Given the primary current density Jp, calculate the electric potential on the surface of the medium. On the contrary, calculating J<sup>p</sup> from the knowledge of U consists of the inverse EEG problem.

So far, we managed to derive at an equation which allows the computation of the electric potential, but our framework is still incomplete. Because the medium under consideration is finite, that is, confined in space by a closed surface-boundary, we need a set of additional constrains, the so-called boundary conditions, which U has to satisfy as well. This stems from the fact that when the medium, in which a wave propagates, displays alterations in its material properties (e.g. different conductivities), the wave is reflected, transmitted or both. In any case, at the interface S separating two regions V1, V<sup>2</sup> the wave and its normal derivative must be continuous at the interface, namely

$$\mathcal{U}\mathcal{U}\_1 = \mathcal{U}\_2 \text{ and } \sigma\_1 \eth\_\nu \mathcal{U}\_1 = \sigma\_2 \eth\_\nu \mathcal{U}\_2 \text{ on } \mathcal{S}, \tag{12}$$

where σ1, σ<sup>2</sup> denote the conductivity of V1, V2, respectively. Above relations are valid only if no charged layer near the interface exists, that is, the absence of primary currents in the vicinity of S is secured. The first of Eq. (12) is known as Dirichlet condition, whereas the second is called a Neumann condition. Neumann conditions must be supplemented by the compatibility condition that the sum of all contributions of the normal derivative ∂νU on S must cancel out, that is,

$$
\oint \partial\_{\nu} \mathcal{U} \, \mathbf{d} \, \mathbf{S} = 0. \tag{13}
$$

Further note that, by virtue of ∂νU, the value of the electric potential is non-unique up to an additive constant. Eqs. (11)–(13) are all we need in order to tackle the EEG problem.

### 3. The brain modelled as a volume conductor

The next step in our journey is to introduce a geometry simple enough in order to carry out the mathematics associated, still adequate realistic in order to illustrate what happens. In practice, above specifications are never met. On the grounds that our interest is focused on the derivation of analytic formulas which will allow us to identify and recognize underlying phenomena, we restrain ourselves to the study of two particular geometries: (i) the spherical and (ii) the ellipsoidal. The distinctness of these two geometries lies in a different representation of the same point in space.

Before we continue, we have to incorporate our assumption regarding the nature of the primary current. For a single, localized dipole at point r<sup>0</sup> and moment Q, we consider

$$\mathbf{J}\_{\mathbf{p}} = \mathbf{Q}\delta(\mathbf{r} - \mathbf{r}\_0),\tag{14}$$

where the functional δðrÞ is the Dirac measure, a concept included in the analysis by the reason of representing the concentration of J<sup>p</sup> to a single point [21].

#### 3.1. The homogeneous spherical brain

on the surface of the medium. On the contrary, calculating J<sup>p</sup> from the knowledge of U consists

So far, we managed to derive at an equation which allows the computation of the electric potential, but our framework is still incomplete. Because the medium under consideration is finite, that is, confined in space by a closed surface-boundary, we need a set of additional constrains, the so-called boundary conditions, which U has to satisfy as well. This stems from the fact that when the medium, in which a wave propagates, displays alterations in its material properties (e.g. different conductivities), the wave is reflected, transmitted or both. In any case, at the interface S separating two regions V1, V<sup>2</sup> the wave and its normal derivative must be

where σ1, σ<sup>2</sup> denote the conductivity of V1, V2, respectively. Above relations are valid only if no charged layer near the interface exists, that is, the absence of primary currents in the vicinity of S is secured. The first of Eq. (12) is known as Dirichlet condition, whereas the second is called a Neumann condition. Neumann conditions must be supplemented by the compatibility condition that the sum of all contributions of the normal derivative ∂νU on S must cancel out,

Further note that, by virtue of ∂νU, the value of the electric potential is non-unique up to an

The next step in our journey is to introduce a geometry simple enough in order to carry out the mathematics associated, still adequate realistic in order to illustrate what happens. In practice, above specifications are never met. On the grounds that our interest is focused on the derivation of analytic formulas which will allow us to identify and recognize underlying phenomena, we restrain ourselves to the study of two particular geometries: (i) the spherical and (ii) the ellipsoidal. The distinctness of these two geometries lies in a different representation of the

Before we continue, we have to incorporate our assumption regarding the nature of the

where the functional δðrÞ is the Dirac measure, a concept included in the analysis by the reason

primary current. For a single, localized dipole at point r<sup>0</sup> and moment Q, we consider

additive constant. Eqs. (11)–(13) are all we need in order to tackle the EEG problem.

3. The brain modelled as a volume conductor

of representing the concentration of J<sup>p</sup> to a single point [21].

U<sup>1</sup> ¼ U<sup>2</sup> and σ1∂νU<sup>1</sup> ¼ σ2∂νU<sup>2</sup> on S, ð12Þ

∮ ∂νU dS ¼ 0: ð13Þ

J<sup>p</sup> ¼ Qδðr � r0Þ, ð14Þ

of the inverse EEG problem.

40 Electroencephalography

continuous at the interface, namely

that is,

same point in space.

The simplest possible geometry in order to represent the brain consists of a spherical homogeneous conductor with radius a and conductivity σ. The reasons are twofold. First of all, it is a good fit to the actual brain. Secondly, it is the only geometry for which data regarding geometrical and physical aspects are immediately available. For example, knowledge of the brains volume can be directly translated into the radius of the corresponding sphere. Importantly, the spherical geometry allows the deduction of explicit expressions for the quantities involved and therefore permits a thorough investigation of the behaviour of the system under scrutiny without running every single time a—mostly time-consuming—computer simulation.

The obvious question at hand is what part of the brain do we model? Clearly, since our initial model is based on homogeneity, it cannot represent the brain en masse, as mentioned in Section 1. So we start with the uppermost region of the human central nervous system, the cerebrum. It is therefore of uttermost importance to be aware of the strengths and weaknesses of the proposed model(s). Without any doubt, the homogeneous spherical model presents an unrealistic assumption of the brain-head system. So why should we bother with theoretical models? The answer is relatively simple. We need them in order to be able to draw conclusions when we move to build models of higher complexity. They serve the important task revealing gaps between forthcoming models, but more substantially they allow us to test the reliability of the introduced algorithms in a straightforward and timely matter. On the other hand, with the homogeneous model, activity in subcortical structures is impossible to detect. Moreover, the influence of the bone architecture enclosing the brain cannot be assessed as well. The latter implies that we actually do not record EEG data at all, but rather monitor the electrophysiology of the (exposed) cerebrum by electrocorticography (ECoG), or intracranial electroencephalography (iEEG).

#### 3.1.1. Forward and inverse problem for a single dipole and multiple dipoles

Having aforementioned remarks in mind, let's begin finding a relation which connects the electric potential on the surface of our conductor model with the electric activity of cells inside. Our goal is achieved solving Eq. (11) combined with expression (14), namely

$$
\Delta \mathcal{U} = \sigma^{-1} \mathbf{Q} \cdot \nabla \delta(\mathbf{r} - \mathbf{r}\_0), \; r < a \tag{15}
$$

supplemented by the condition

$$
\partial\_r \mathcal{U} = 0, \; r = a,\tag{16}
$$

which follows at once from the second condition (12) expressing the circumstance that the conductivity outside the brain vanishes, and expresses the 'reality' that no current exists outside the brain.<sup>1</sup> Note that the compatibility condition (13) is automatically satisfied by Eq. (16).

<sup>1</sup> This statement is true only for a homogeneous conductor as described in Section 3.

Employing analytic techniques, the interested reader can find all details in Ref. [22]; it is not hard to show that the solution regarding Eqs. (15) and (16) evaluated at the surface is given as

$$\mathbf{M}\_{\text{surf}}(\hat{\mathbf{r}}, \mathbf{r}\_0) = \sum\_{n=1}^{\infty} \sum\_{m=-n}^{n} \mathbf{A}\_n^{\text{m}}(\hat{\mathbf{r}}\_0) Y\_n^{\text{m}}(\hat{\mathbf{r}}), \ \mathbf{A}\_n^{\text{m}}(\hat{\mathbf{r}}\_0) = \frac{\mathbf{Q}}{\sigma n a^{n+1}} \cdot \nabla\_{\mathbf{r}\_0} \left( r\_0^n \overline{Y}\_n^{\text{m}}(\hat{\mathbf{r}}\_0) \right), \tag{17}$$

where a hat '^'denotes the unit vector, Y and Y are the spherical harmonics and corresponding complex conjugate, respectively, which are complex-valued functions, the latter with equal real part and imaginary part equal in magnitude but opposite in sign [23].

Expression (17) can be simplified using a summation formula [22, 24] yielding the following closed form:

$$\mathcal{U}\_{\rm surf}(\hat{\mathbf{r}}, \mathbf{r}\_0) = \frac{\mathbf{Q}}{4\pi\sigma} \cdot \left(2\frac{\mathbf{R}}{R^3} + \frac{1}{aR}\frac{R\hat{\mathbf{r}} + \mathbf{R}}{R + \hat{\mathbf{r}} \cdot \mathbf{R}}\right), \; \mathbf{R} = a\hat{\mathbf{r}} - \mathbf{r}\_0. \tag{18}$$

where an italics-type capital letter denotes the magnitude of the corresponding vector. Eq. (18) consists of the simplest, straightforward expression regarding EEG data.

Let us now concentrate on the most important aspect when it comes to imaging modalities, such as EEG, namely the problem of identifying the primary source by means of a generated electromagnetic field. We recall that the notion of an equivalent dipole source has been adopted in order to summarize the entire microscopic currents located in the vicinity of a specific area in the brain. Notwithstanding, there does not exist an exclusive source configuration for each set of electroencephalographic measurements, constituting the corresponding inverse problem non-unique. The only way to eliminate non-uniqueness is to provide supplementary information's, that is, imposing additional assumptions. By introducing the assumption of an equivalent dipole source, the inverse problem can be solved exactly as we will show in the sequel.

The inverse problem for the homogeneous spherical conductor is formulated as follows. From surface measurements, we identify the potential Usurf, given via Eq. (17), from which we have to calculate the position r<sup>0</sup> ¼ ðx0, y0, z0Þ and moment Q ¼ ðQx, Qy, QzÞ of the dipole. In our case, this is easily achieved, noting that every detail regarding r<sup>0</sup> and Q are encoded into the coefficients A<sup>m</sup> <sup>n</sup> of expansion (17). The latter are evaluated with the aid of the orthogonality condition

$$\mathbf{A}\_{n}^{\prime\prime} = \oint I\_{\text{surf}}(\hat{\mathbf{r}} \,\prime \mathbf{r}\_{0}) \overline{Y}\_{n}^{\prime\prime}(\hat{\mathbf{r}}) \,\mathrm{d}S(\hat{\mathbf{r}}).\tag{19}$$

A word of caution with respect to Eq. (19). In order to employ the latter, we must know the surface potential on every single point via Eq. (17). In practice, the function Usurf is acquired as a continuous function via interpolation of the discrete set of EEG measurements. Moreover, since each vector contains three coordinates, in order to pinpoint the position and moment of the dipole, at least six equations are required. As a result, we expand the coefficients A<sup>m</sup> <sup>n</sup> , given via Eq. (17), for n ¼ 1 and 2 providing eight relations in total. The first three of them, namely A�<sup>1</sup> <sup>1</sup> , <sup>A</sup><sup>0</sup> <sup>1</sup> and A<sup>1</sup> <sup>1</sup>, are proportional to the components Qx, Qy, Qz of Q, and thus we find

Employing analytic techniques, the interested reader can find all details in Ref. [22]; it is not hard to show that the solution regarding Eqs. (15) and (16) evaluated at the surface is

<sup>n</sup> <sup>ð</sup>^rÞ, <sup>A</sup><sup>m</sup>

where a hat '^'denotes the unit vector, Y and Y are the spherical harmonics and corresponding complex conjugate, respectively, which are complex-valued functions, the latter with equal real

Expression (17) can be simplified using a summation formula [22, 24] yielding the following

1 aR

where an italics-type capital letter denotes the magnitude of the corresponding vector. Eq. (18)

Let us now concentrate on the most important aspect when it comes to imaging modalities, such as EEG, namely the problem of identifying the primary source by means of a generated electromagnetic field. We recall that the notion of an equivalent dipole source has been adopted in order to summarize the entire microscopic currents located in the vicinity of a specific area in the brain. Notwithstanding, there does not exist an exclusive source configuration for each set of electroencephalographic measurements, constituting the corresponding inverse problem non-unique. The only way to eliminate non-uniqueness is to provide supplementary information's, that is, imposing additional assumptions. By introducing the assumption of an equivalent dipole source, the inverse problem can be solved exactly as we will show

The inverse problem for the homogeneous spherical conductor is formulated as follows. From surface measurements, we identify the potential Usurf, given via Eq. (17), from which we have to calculate the position r<sup>0</sup> ¼ ðx0, y0, z0Þ and moment Q ¼ ðQx, Qy, QzÞ of the dipole. In our case, this is easily achieved, noting that every detail regarding r<sup>0</sup> and Q are encoded into the

<sup>n</sup> <sup>¼</sup> <sup>∮</sup> <sup>U</sup>surfð^r,r0ÞY<sup>m</sup>

A word of caution with respect to Eq. (19). In order to employ the latter, we must know the surface potential on every single point via Eq. (17). In practice, the function Usurf is acquired as a continuous function via interpolation of the discrete set of EEG measurements. Moreover, since each vector contains three coordinates, in order to pinpoint the position and moment of the dipole, at least six equations are required. As a result, we expand the coefficients A<sup>m</sup>

<sup>n</sup> of expansion (17). The latter are evaluated with the aid of the orthogonality

� �

R^r þ R R þ ^r � R

<sup>R</sup><sup>3</sup> <sup>þ</sup>

<sup>n</sup> <sup>ð</sup>^r0Þ ¼ <sup>Q</sup>

<sup>σ</sup>nanþ<sup>1</sup> � <sup>∇</sup><sup>r</sup><sup>0</sup>

� r n 0Y<sup>m</sup> <sup>n</sup> ð^r0Þ �

, R ¼ a^r � r0, ð18Þ

<sup>n</sup> ð^rÞdSð^rÞ: ð19Þ

<sup>n</sup> , given

, ð17Þ

given as

42 Electroencephalography

closed form:

in the sequel.

coefficients A<sup>m</sup>

condition

<sup>U</sup>surfð^r,r0Þ ¼ <sup>X</sup><sup>∞</sup>

n¼1

<sup>U</sup>surfð^r,r0Þ ¼ <sup>Q</sup>

Xn m¼�n

A<sup>m</sup> <sup>n</sup> <sup>ð</sup>^r0ÞY<sup>m</sup>

part and imaginary part equal in magnitude but opposite in sign [23].

<sup>4</sup>πσ � <sup>2</sup> <sup>R</sup>

consists of the simplest, straightforward expression regarding EEG data.

A<sup>m</sup>

$$Q\_x = \sqrt{\frac{2\pi}{3}} \sigma a^2 (\mathbf{A}\_1^{-1} - \mathbf{A}\_1^1), \ Q\_y = -i\sqrt{\frac{2\pi}{3}} \sigma a^2 (\mathbf{A}\_1^{-1} + \mathbf{A}\_1^1), \ Q\_z = 2\sqrt{\frac{\pi}{3}} \sigma a^2 \mathbf{A}\_1^0 \tag{20}$$

where i denotes the imaginary unit. It is not hard to show that whereas the difference A�<sup>1</sup> <sup>1</sup> � <sup>A</sup><sup>1</sup> 1 is real valued, the corresponding sum is imaginary. The remaining five coefficients A<sup>m</sup> <sup>2</sup> for n ¼ 2 and m ¼ �2,�1,0,1,2 connect the position of the dipole with its moments which are eliminated utilizing Eq. (20). Again, after some algebra we have

$$x\_0 = \frac{a}{\sqrt{5}} \left( \frac{\mathbf{A}\_2^{-2}}{\mathbf{A}\_1^{-1}} - \frac{\mathbf{A}\_2^2}{\mathbf{A}\_1^1} \right), \ y\_0 = -i \frac{a}{\sqrt{5}} \left( \frac{\mathbf{A}\_2^{-2}}{\mathbf{A}\_1^{-1}} + \frac{\mathbf{A}\_2^2}{\mathbf{A}\_1^1} \right), \ z\_0 = \frac{2a}{\sqrt{5}} \frac{1}{\mathbf{A}\_1^1} \left( \mathbf{A}\_2^1 - \frac{\mathbf{A}\_2^2 \mathbf{A}\_1^0}{\sqrt{2} \mathbf{A}\_1^1} \right). \tag{21}$$

There is another way to recover the solution to the inverse problem for a single dipole. The approach illustrated provides a glimpse into the beauty of mathematical analysis. We will show that the uniqueness of the inverse problem, for a single dipole, is closely connected with the condition of attaining certain relations connecting the measured data. Considering that we need six equations to identify the source, we expand Eq. (17) for n ¼ 1, 2 and express the resulting relation in Cartesian coordinates, yielding [25]

$$\mathbf{M}\_{\text{surf}} = A\_1 \mathbf{x}\_0 + A\_2 \mathbf{y}\_0 + A\_3 \mathbf{z}\_0 + B\_1 \mathbf{x}\_0^2 + B\_2 \mathbf{y}\_0^2 + B\_3 \mathbf{z}\_0^2 + C\_1 \mathbf{x}\_0 \mathbf{y}\_0 + C\_2 \mathbf{x}\_0 \mathbf{z}\_0 + C\_3 \mathbf{y}\_0 \mathbf{z}\_0 + \dots \quad (22)$$

where again the first three coefficients A1, A2, A<sup>3</sup> are proportional to the components of Q, whereas the remaining six involve products of the position of the dipole with its moments. Substituting the expressions for Qx, Qy, Qz into B1, B2, B3, C1, C2, C3, we arrive at the linear but overdetermined system Ar<sup>0</sup> ¼ b with

$$\mathbf{A} = \begin{bmatrix} 2A\_1 & -A\_2 & -A\_3 \\ -A\_1 & 2A\_2 & -A\_3 \\ -A\_1 & -A\_2 & 2A\_3 \\ A\_2 & A\_1 & 0 \\ A\_3 & 0 & A\_1 \\ 0 & A\_3 & A\_2 \end{bmatrix}, \mathbf{b} = \frac{6a^2}{5} \begin{bmatrix} B\_1 \\ B\_2 \\ B\_3 \\ C\_1 \\ C\_2 \\ C\_3 \end{bmatrix}. \tag{23}$$

Adopting Gauss elimination, we find that the position of the dipole is given as

$$\mathbf{x}\_{0} = \frac{3}{2A\_{1}A\_{2}A\_{3}}(A\_{1}\mathbf{a}) \cdot \tilde{\mathbf{b}},\\\ y\_{0} = \frac{3}{2A\_{1}A\_{2}A\_{3}}(A\_{2}\mathbf{a}) \cdot (\mathbb{R}\_{z}\tilde{\mathbf{b}}),\\\ y\_{0} = \frac{3}{2A\_{1}A\_{2}A\_{3}}(A\_{3}\mathbf{a}) \cdot (\mathbb{R}\_{y}\tilde{\mathbf{b}}),\qquad(24)$$

where a ¼ ðA1,A2,A3Þ T, <sup>b</sup><sup>~</sup> ¼ ð�C1,C2,C3<sup>Þ</sup> T, T denotes transposition, whereas ℝ<sup>y</sup> and ℝ<sup>z</sup> denote the rotation matrices about the y- and z-axis, respectively.

A very important aspect, revealed by the above analysis, is evidence that the dipoles position specified by relations (24) is unique, only if the recorded values for the coefficients present in Eq. (22) satisfy the following relations:

$$\begin{aligned} 2A\_1A\_2A\_3(B\_3-B\_2) + 3A\_1(A\_2^2-A\_3^2)\mathbb{C}\_2 + 3(A\_2^2+A\_3^2)(A\_3\mathbb{C}\_1-A\_2\mathbb{C}\_3) &= 0, \\ -2A\_1A\_2A\_3(B\_2+2B\_3) + 3A\_2(A\_3^2-A\_1^2)\mathbb{C}\_3 + 3(A\_1^2+A\_3^2)(A\_1\mathbb{C}\_2-A\_3\mathbb{C}\_1) &= 0. \end{aligned} \tag{25}$$

Indeed, replacing the analytic expressions for the coefficients of Eq. (22) into Eq. (25), these are trivially satisfied. Complementary, we briefly state that the least-square solution to Eq. (23) is r<sup>0</sup> ¼ A<sup>þ</sup>b þ ðI � A<sup>þ</sup>AÞY, where A<sup>þ</sup> is the pseudoinverse of A, and Y is an arbitrary vector [26]. Evidently, the presence of Y constitutes the solution non-unique. Nonetheless, for a single dipole the equality A<sup>þ</sup>A ¼ I holds true, leading to the minimum norm solution r<sup>0</sup> ¼ A<sup>þ</sup>b:

However, in 'reality' EEG recordings provide values for the coefficients A<sup>m</sup> <sup>n</sup> , as given in Eq. (17), but our uniqueness criteria (25) are bound to the coefficients of Eq. (22). So, how do the correlate? Well, connection formulas are established by expanding the first eight terms of the surface potential (17). These terms are then rearranged to form an expression similar to Eq. (22). Comparing the resulting relation with Eq. (22) provides the connection. The final formulas can be found in Ref. [25].

Until now, we attended the situation where brain activity is simulated by a single dipole. What happens if a larger area or multiple areas are stimulated? Is a single dipole adequate to describe the event? Multiple areas of the cortex are often expected to be active at the same time, so the answer must be no. In Section 2, we mentioned that the electric fields are superposable. As a consequence, the surface potential due to N-dipole sources is computed as the sum of the potentials Uj corresponding to dipoles located at r0<sup>j</sup> and strength Qj: In order to calculate their positions and moments, we require at least 6N equations. Although an analytic inversion algorithm can be derived [27], the steps necessary are cumbersome and beyond the scope of the present document. It remains, however, the question of identifying multiple localized sources by means of a recorded potential. Is it possible to be led into an erroneous conclusion when we have to recognize the number of activated areas? This situation can occur when data (coefficients) received are falsely interpreted as evoked by a single dipole. To show this, we first rewrite the surface potential resulting from a single dipole (17), in the form

$$\mathcal{U}\_{\rm surf}(\hat{\mathbf{r}}, \mathbf{r}\_0) = \frac{1}{4\pi\sigma} \sum\_{n=1}^{\ast} \frac{2n+1}{na^{n+1}} (\mathbf{Q} \cdot \nabla\_{\mathbf{r}\_0}) r\_0^n P\_n(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}\_0), \tag{26}$$

employing the addition theorem for Legendre functions Pn: The same expression is valid on the occasion of N dipoles, that is,

$$\mathcal{U}\_{\text{surf}}(\hat{\mathbf{r}}, \mathbf{r}\_{0\uparrow}) = \frac{1}{4\pi\sigma} \sum\_{n=1}^{\infty} \frac{2n+1}{n a^{n+1}} \sum\_{j=1}^{N} (\mathbf{Q}\_j \cdot \nabla\_{\mathbf{r}\_{0\uparrow}}) r\_{0\uparrow}^n P\_n(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}\_{0\uparrow}). \tag{27}$$

Since we consider the surface potentials to be identical, irrespective of the number of dipoles, the coefficients of Eqs. (26) and (27) must be equal, namely

$$p(\mathbf{Q} \cdot \nabla\_{\mathbf{r}\_0}) r\_0^\mu P\_n(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}\_0) = \sum\_{j=1}^N (\mathbf{Q}\_j \cdot \nabla\_{\mathbf{r}\_0}) r\_{0j}^\mu P\_n(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}\_{0j}).\tag{28}$$

If Eq. (28) is satisfied, a distinguishing between dipoles is not possible. A thorough investigation of the latter shows that if the dipoles are parallel to each other, the uniqueness conditions (25) are fulfilled and it is impossible to decide if the measurements are induced by a single or finite number of dipoles [25].

#### 3.1.2. Forward and inverse problem for distributed activity

<sup>2</sup>A1A2A3ðB<sup>3</sup> � <sup>B</sup>2Þ þ <sup>3</sup>A1ðA<sup>2</sup>

44 Electroencephalography

formulas can be found in Ref. [25].

the occasion of N dipoles, that is,

�2A1A2A3ðB<sup>2</sup> <sup>þ</sup> <sup>2</sup>B3Þ þ <sup>3</sup>A2ðA<sup>2</sup>

<sup>2</sup> � <sup>A</sup><sup>2</sup>

However, in 'reality' EEG recordings provide values for the coefficients A<sup>m</sup>

<sup>U</sup>surfð^r,r0Þ ¼ <sup>1</sup>

<sup>U</sup>surfð^r,r0jÞ ¼ <sup>1</sup>

the coefficients of Eqs. (26) and (27) must be equal, namely

ðQ � ∇<sup>r</sup><sup>0</sup> Þr

n

4πσ

4πσ

X∞ n¼1

<sup>0</sup>Pnð^<sup>r</sup> � ^r0Þ ¼ <sup>X</sup>

X∞ n¼1

2n þ 1

employing the addition theorem for Legendre functions Pn: The same expression is valid on

2n þ 1 nanþ<sup>1</sup>

Since we consider the surface potentials to be identical, irrespective of the number of dipoles,

N

j¼1

X N

ðQ<sup>j</sup> � ∇<sup>r</sup>0<sup>j</sup>

j¼1

ðQ<sup>j</sup> � ∇<sup>r</sup>0<sup>j</sup>

Þr n 0j

nanþ<sup>1</sup> <sup>ð</sup><sup>Q</sup> � <sup>∇</sup><sup>r</sup><sup>0</sup> <sup>Þ</sup><sup>r</sup>

n

Þr n 0j

<sup>0</sup>Pnð^r � ^r0Þ, ð26Þ

Pnð^r � ^r0jÞ: ð27Þ

Pnð^r � ^r0jÞ: ð28Þ

<sup>3</sup> � <sup>A</sup><sup>2</sup>

<sup>3</sup>ÞC<sup>2</sup> <sup>þ</sup> <sup>3</sup>ðA<sup>2</sup>

Indeed, replacing the analytic expressions for the coefficients of Eq. (22) into Eq. (25), these are trivially satisfied. Complementary, we briefly state that the least-square solution to Eq. (23) is r<sup>0</sup> ¼ A<sup>þ</sup>b þ ðI � A<sup>þ</sup>AÞY, where A<sup>þ</sup> is the pseudoinverse of A, and Y is an arbitrary vector [26]. Evidently, the presence of Y constitutes the solution non-unique. Nonetheless, for a single dipole the equality A<sup>þ</sup>A ¼ I holds true, leading to the minimum norm solution r<sup>0</sup> ¼ A<sup>þ</sup>b:

Eq. (17), but our uniqueness criteria (25) are bound to the coefficients of Eq. (22). So, how do the correlate? Well, connection formulas are established by expanding the first eight terms of the surface potential (17). These terms are then rearranged to form an expression similar to Eq. (22). Comparing the resulting relation with Eq. (22) provides the connection. The final

Until now, we attended the situation where brain activity is simulated by a single dipole. What happens if a larger area or multiple areas are stimulated? Is a single dipole adequate to describe the event? Multiple areas of the cortex are often expected to be active at the same time, so the answer must be no. In Section 2, we mentioned that the electric fields are superposable. As a consequence, the surface potential due to N-dipole sources is computed as the sum of the potentials Uj corresponding to dipoles located at r0<sup>j</sup> and strength Qj: In order to calculate their positions and moments, we require at least 6N equations. Although an analytic inversion algorithm can be derived [27], the steps necessary are cumbersome and beyond the scope of the present document. It remains, however, the question of identifying multiple localized sources by means of a recorded potential. Is it possible to be led into an erroneous conclusion when we have to recognize the number of activated areas? This situation can occur when data (coefficients) received are falsely interpreted as evoked by a single dipole. To show this, we first rewrite the surface potential resulting from a single dipole (17), in the form

<sup>1</sup>ÞC<sup>3</sup> <sup>þ</sup> <sup>3</sup>ðA<sup>2</sup>

<sup>2</sup> <sup>þ</sup> <sup>A</sup><sup>2</sup>

<sup>1</sup> <sup>þ</sup> <sup>A</sup><sup>2</sup>

<sup>3</sup>ÞðA3C<sup>1</sup> � A2C3Þ ¼ 0,

<sup>3</sup>ÞðA1C<sup>2</sup> � A3C1Þ ¼ 0:

ð25Þ

<sup>n</sup> , as given in

If we abandon the assumption that the primary current J<sup>p</sup> is represented by a dipole, the additional information leading to a unique solution regarding the inverse problem is automatically lost. Notwithstanding, dropping the dipole hypothesis allows us to simulate complicated activation patterns in terms of distributed currents. The potential on the surface is computed solving Eq. (11) directly, accompanied by proper boundary conditions. For the inverse source problem, non-uniqueness remains a point of concern.

Albanese and Monk [28] illustrated that it is not possible to recreate a three-dimensional current based on EEG measurements. This result has been practically demonstrated in Ref. [27], where the authors show that the radius of a small spherical current cannot be recovered. For currents having dimensions less than three, the inverse EEG problem admits a unique solution. In previous sections, we swiftly examined currents of zero dimensionality, namely dipoles. In what follows, we will explore the forward and inverse EEG problems for one- and two-dimensional continuously distributed currents [29, 30].

We begin by assuming that the current is a small line segment of length 2L, centred at r<sup>0</sup> and oriented along an arbitrary direction α^: The primary current is then approximated as J<sup>p</sup> ≃ Q þ tA, where t is a variable taking values in the interval [-L,L], whereas A is the directional derivative of J<sup>p</sup> along α^: Replacing the approximation of J<sup>p</sup> into Eq. (11) provides the surface potential in the case of a linearly distributed current [29]. Knowledge of the surface measurements enables us to identify the position, moment, orientation and size of the current. Consequently, we have to determine 13 parameters and require a sufficient number of equations to be able to perform the identification. Operationally, the procedure in order to obtain this set of equations is in principle the same as for a single dipole. We expand the surface potential in a series of harmonic, homogeneous polynomials in Cartesian coordinates, where the coefficients of each monomial are known. Moreover, each coefficient contains a certain number of unknown parameters. Here, we must determine 13 unknowns which means we have to analytically calculate at least 13 coefficients, building the necessary system of equations. However, as the number of unknowns grows the concluding system of equations turns highly nonlinear. In the one-dimensional case, the system consists of a total of 19 equations with four constrains. It is possible to solve the latter at least semi-analytically [29].

When investigating the two-dimensional case, the mathematical complexity takes it up a notch. Assuming that the current is a small disk of radius ε, centred at r<sup>0</sup> and perpendicularly oriented to the vector <sup>r</sup>0, the primary current can be approximated as <sup>J</sup><sup>p</sup> <sup>≃</sup> <sup>Q</sup> <sup>þ</sup> <sup>r</sup> � <sup>D</sup><sup>~</sup> : The quantity <sup>D</sup><sup>~</sup> is called dyadic, a second-order tensor containing physical or geometrical information. The surface potential is provided replacing J<sup>p</sup> into Eq. (11) and solving the corresponding boundary value problem. Reaching that goal facilitates long and tedious manipulations involving integration as well [30]. The inverse problem, on the other hand, follows the guidelines outlined earlier. Due to the specific orientation of the primary current, only seven parameters have to be determined. Once more, the system of equations from which these parameters will be decided is nonlinear [30].

An independent view to the particular problem has been provided by Fokas [31]. After formulating the surface potential with explicit Q dependence (26), computing the surface potential for a continuously distributed current is straightforward. One has to replace Q by J<sup>p</sup> and integrate the resulting expression with respect to the volume of the conductor. The associated manipulations can be simplified by introducing Helmholtz decomposition for the primary current, which states that any three-dimensional smooth vector field can be resolved into the sum of an irrotational and solenoidal vector field, ∇Ψðr0Þ and Aðr0Þ, respectively. As a result, assuming J<sup>p</sup> ¼ ∇Ψ þ ∇ � A under the condition that ∇ � A ¼ 0, we find

$$dU\_{\rm surf}(\hat{\mathbf{r}}, \mathbf{r}\_0) = \oint \sum\_{n=1}^{\infty} \mathbf{C}\_n^m(\Delta \Psi) r\_0^n P\_n(\hat{\mathbf{r}} \cdot \hat{\mathbf{r}}\_0) \mathbf{d}v(\hat{\mathbf{r}}\_0). \tag{29}$$

Note that the electric potential depends only on ΔΨðr0Þ. The integral in Eq. (29), obtained using Gauss theorem and integration by parts, can be computed analytically if we expand Ψðr0Þ in terms of spherical harmonics, namely

$$\Psi(\mathbf{r}\_0) = \sum\_{n=1}^{\infty} \sum\_{m=-n}^{n} \psi\_n^m(r\_0) Y\_n^m(\mathbf{\hat{r}}\_0),\tag{30}$$

furnishing

$$\delta L\_{\text{surf}} = \sum\_{n=1}^{\infty} \sum\_{m=-n}^{n} (\alpha\_n \dot{\psi}\_n^m - \beta\_n \psi\_n^m) Y\_n^m(\hat{\mathbf{r}}). \tag{31}$$

It is remarkable that the above coefficients can be determined only under the assumption that Ψðr0Þ is a bi-harmonic function, namely a solution to the bi-harmonic operator ΔΔ: However, the precise description of the coefficients ψ<sup>m</sup> <sup>n</sup> remains open.

#### 3.2. The homogeneous ellipsoidal brain

From the point of view of mathematical analysis, any three-dimensional object, such as the brain, would be best approximated with the aid of a coordinate system with three degrees of freedom, one in each direction. Fortunately, aforesaid system exists and is called the ellipsoidal coordinate system. Whereas the spherical coordinate system consists of concentric spheres centred at the origin, in the ellipsoidal coordinate system, each point is specified by the intersection of three non-degenerate second-degree surfaces, corresponding to an ellipsoid, a hyperboloid of one sheet as well as a hyperboloid of two sheets. This constitutes the ellipsoidal system significantly more complex and demanding than the spherical one (see [32] for an analytic account). For example, whereas knowledge of the brain's volume can be directly translated into the radius of the corresponding sphere, the volume of an ellipsoid is proportional to the product of three parameters a1, a2, a3, the so-called semi-axes of the ellipsoid. Hence, there does not exist a unique combination of those three parameters providing the volume of the brain. In spite of the

ramifications, the ellipsoidal system is an environment, which allows the installation and interpretation of analytical algorithms to a great extent as well. This stems from the fact that it is the most general system where the Laplacian operator assumes a spectral decomposition. As we dive into this particular system, we recognize that the ellipsoidal geometry is responsible for drastic variations in the behaviour of EEG when compared to the sphere.

#### 3.2.1. Forward and inverse problem for a single dipole and distributed activity

The inverse problem, on the other hand, follows the guidelines outlined earlier. Due to the specific orientation of the primary current, only seven parameters have to be determined. Once more, the

An independent view to the particular problem has been provided by Fokas [31]. After formulating the surface potential with explicit Q dependence (26), computing the surface potential for a continuously distributed current is straightforward. One has to replace Q by J<sup>p</sup> and integrate the resulting expression with respect to the volume of the conductor. The associated manipulations can be simplified by introducing Helmholtz decomposition for the primary current, which states that any three-dimensional smooth vector field can be resolved into the sum of an irrotational and solenoidal vector field, ∇Ψðr0Þ and Aðr0Þ, respectively. As a

system of equations from which these parameters will be decided is nonlinear [30].

result, assuming J<sup>p</sup> ¼ ∇Ψ þ ∇ � A under the condition that ∇ � A ¼ 0, we find

<sup>Ψ</sup>ðr0Þ ¼ <sup>X</sup><sup>∞</sup>

n¼1

<sup>U</sup>surf <sup>¼</sup> <sup>X</sup><sup>∞</sup>

n¼1

Xn m¼�n

X∞ n¼1 Cm <sup>n</sup> ðΔΨÞr n

Note that the electric potential depends only on ΔΨðr0Þ. The integral in Eq. (29), obtained using Gauss theorem and integration by parts, can be computed analytically if we expand Ψðr0Þ in

> Xn m¼�n

> > <sup>ð</sup>αnψ\_ <sup>m</sup>

It is remarkable that the above coefficients can be determined only under the assumption that Ψðr0Þ is a bi-harmonic function, namely a solution to the bi-harmonic operator ΔΔ: However,

From the point of view of mathematical analysis, any three-dimensional object, such as the brain, would be best approximated with the aid of a coordinate system with three degrees of freedom, one in each direction. Fortunately, aforesaid system exists and is called the ellipsoidal coordinate system. Whereas the spherical coordinate system consists of concentric spheres centred at the origin, in the ellipsoidal coordinate system, each point is specified by the intersection of three non-degenerate second-degree surfaces, corresponding to an ellipsoid, a hyperboloid of one sheet as well as a hyperboloid of two sheets. This constitutes the ellipsoidal system significantly more complex and demanding than the spherical one (see [32] for an analytic account). For example, whereas knowledge of the brain's volume can be directly translated into the radius of the corresponding sphere, the volume of an ellipsoid is proportional to the product of three parameters a1, a2, a3, the so-called semi-axes of the ellipsoid. Hence, there does not exist a unique combination of those three parameters providing the volume of the brain. In spite of the

<sup>n</sup> remains open.

<sup>n</sup> � <sup>β</sup>nψ<sup>m</sup>

<sup>n</sup> <sup>Þ</sup>Y<sup>m</sup>

ψ<sup>m</sup> <sup>n</sup> <sup>ð</sup>r0ÞY<sup>m</sup>

<sup>0</sup>Pnð^r � ^r0Þdvð^r0Þ: ð29Þ

<sup>n</sup> ð^r0Þ, ð30Þ

<sup>n</sup> ð^rÞ: ð31Þ

<sup>U</sup>surfð^r,r0Þ ¼ ∮

terms of spherical harmonics, namely

the precise description of the coefficients ψ<sup>m</sup>

3.2. The homogeneous ellipsoidal brain

furnishing

46 Electroencephalography

In order to acquire the surface potential for an active dipole within an ellipsoidal brain, the same machinery as for the sphere is utilized. There are, however, some differences. The spherical coordinates are easily established by fixing the centre of the system and moving a distance r away. Symmetry will do the rest. In the ellipsoidal system, the procedure works a lot different. In order to solve boundary value problems in the ellipsoidal coordinate system, such as the forward EEG problem, it is essential to adopt an ellipsoid in such a way as to fit the actual boundary of the conductor under consideration, by choosing a particular value of the ellipsoidal 'radial' variable ρ. This is secured if we use the boundary of our domain to be the socalled reference ellipsoid and construct the ellipsoidal system that is based on it. Further, denote by h1, h2, h<sup>3</sup> the three semifocal distances from which the orthogonal ellipsoidal coordinate system ðρ, ν, μÞ is derived. Each confocal ellipsoidal surface is defined by a constant value of the 'radial' variable ρ∈ ðh2, ∞Þ, with ρ ¼ a<sup>1</sup> indicating the reference ellipsoid.

Consider in what follows a homogeneous ellipsoidal conductor with semi-axes ai, i ¼ 1,2,3 and conductivity σ which takes the place of the reference ellipsoid for our ellipsoidal coordinate system. The surface potential is again computed by solving Eq. (11) combined with expression (14), namely

$$
\Delta \mathcal{U} = \sigma^{-1} \mathbf{Q} \cdot \nabla \delta(\mathbf{r} - \mathbf{r}\_0)\_\prime \ h\_2 < \rho < a\_1. \tag{32}
$$

$$
\partial\_{\rho} \mathcal{U} = 0, \ \rho = a\_1. \tag{33}
$$

Employing analytic techniques, the interested reader can find all details in Ref. [33], the solution regarding Eqs. (29) and (30) evaluated at the surface is

$$\mathcal{U}\_{\text{surf}}(\mathbf{r};\mathbf{r}\_{0}) = \frac{1}{\sigma} \sum\_{n=1}^{\approx} \sum\_{m=1}^{2n+1} \mathbf{B}\_{n}^{m}(\mathbf{r}\_{0}) \mathbf{E}\_{n}^{m}(\mu) \mathbf{E}\_{n}^{m}(\nu), \ \mathbf{B}\_{n}^{m}(\mathbf{r}\_{0}) = \frac{\mathbf{Q} \cdot \nabla\_{\mathbf{r}\_{0}} \mathbf{E}\_{n}^{m}(\mathbf{r}\_{0})}{d}. \tag{34}$$

where E<sup>m</sup> <sup>n</sup> are the Lamé functions, E<sup>m</sup> <sup>n</sup> ðrÞ symbolizes the triple product of these Lamé functions, and d some constant. Compared to the corresponding solutions for the spherical conductor, the above formula looks very similar to Eq. (17). As usual, the devil is in the detail.

Before proceeding, we mention that an elegant and straightforward expression connecting the surface potential with the moment Q and position r<sup>0</sup> of the dipole, similar to Eq. (18), does not exist for the ellipsoidal conductor.

We turn now to the inverse problem. In principle, the procedure described for the sphere is generally applicable since it is geometry independent. Note again that every detail regarding the moment Q and position r<sup>0</sup> of the dipole is encoded into the coefficients of expansion (31). As a result, we expand the coefficients B<sup>m</sup> <sup>n</sup> for n ¼ 1 and 2 yielding eight relations in total. Defining a new parameter gm <sup>n</sup> as

$$\mathbf{g}\_n^m = \mathbf{Q} \cdot \nabla\_{\mathbf{r}\_0} \mathbb{E}\_n^m(\mathbf{r}\_0) = d \, \mathbf{B}\_n^m \tag{35}$$

which incorporates characteristics, such as the geometrical and physical properties of the conductor, as well as the EEG measurements, the algebraic manipulations leading to the solution are somewhat a little more painless. Solving a linear system of six equations with six unknowns, it can be shown [34] that the dipoles position r<sup>0</sup> and moment Q depends only on g1 1, g<sup>2</sup> 1, g<sup>3</sup> 1, g<sup>3</sup> 2, g<sup>4</sup> 2, g<sup>5</sup> <sup>2</sup> as

$$\mathbf{x}\_{0} = \frac{1}{2h\_{2}h\_{3}} \left( \frac{\mathbf{g}\_{2}^{3}}{\mathbf{g}\_{1}^{2}} + \frac{\mathbf{g}\_{2}^{4}}{\mathbf{g}\_{1}^{3}} - \frac{\mathbf{g}\_{1}^{1}\mathbf{g}\_{2}^{5}}{\mathbf{g}\_{1}^{2}\mathbf{g}\_{1}^{3}} \right), \ y\_{0} = \frac{1}{2h\_{1}h\_{3}} \left( \frac{\mathbf{g}\_{2}^{3}}{\mathbf{g}\_{1}^{1}} + \frac{\mathbf{g}\_{2}^{5}}{\mathbf{g}\_{1}^{3}} - \frac{\mathbf{g}\_{1}^{2}\mathbf{g}\_{2}^{4}}{\mathbf{g}\_{1}^{1}\mathbf{g}\_{1}^{3}} \right), \ z\_{0} = \frac{1}{2h\_{1}h\_{2}} \left( \frac{\mathbf{g}\_{2}^{4}}{\mathbf{g}\_{1}^{1}} + \frac{\mathbf{g}\_{2}^{5}}{\mathbf{g}\_{1}^{2}} - \frac{\mathbf{g}\_{1}^{3}\mathbf{g}\_{2}^{3}}{\mathbf{g}\_{1}^{1}\mathbf{g}\_{1}^{2}} \right) \tag{36}$$

and

$$Q\_m = \frac{h\_m}{h\_1 h\_2 h\_3} g\_{1'}^{m} \tag{37}$$

for every <sup>m</sup> <sup>¼</sup> 1,2,3, respectively. Expressing two of the triple products <sup>E</sup><sup>m</sup> <sup>n</sup> ðrÞ in Cartesian coordinates, namely for n ¼ 2 and m ¼ 1,2, we obtain two constrains similar to Eq. (25), which are satisfied by the constants g<sup>1</sup> 1, g<sup>2</sup> 1, g<sup>3</sup> 1, g<sup>3</sup> 2, g<sup>4</sup> 2, g<sup>5</sup> <sup>2</sup> as

$$\mathbf{g}\_2^1 = 2(\Lambda - a\_1^2)(\Lambda - a\_2^2)(\Lambda - a\_3^2) \left( \frac{Q\_1 \mathbf{x}\_0}{\Lambda - a\_1^2} + \frac{Q\_2 \mathbf{y}\_0}{\Lambda - a\_2^2} + \frac{Q\_3 \mathbf{z}\_0}{\Lambda - a\_3^2} \right) \tag{38}$$

and

$$\mathbf{g}\_2^2 = 2(\Lambda' - a\_1^2)(\Lambda' - a\_2^2)(\Lambda' - a\_3^2)\left(\frac{\mathbf{Q}\_1 \mathbf{x}\_0}{\Lambda' - a\_1^2} + \frac{\mathbf{Q}\_2 \mathbf{y}\_0}{\Lambda' - a\_2^2} + \frac{\mathbf{Q}\_3 \mathbf{z}\_0}{\Lambda' - a\_3^2}\right),\tag{39}$$

where Λ, Λ<sup>0</sup> are ellipsoidal parameters depending on the semi-axes ai, i ¼ 1,2,3.

In the case of a continuously distributed current, we follow the example laid out for the sphere. Replacing Q by the primary current J<sup>p</sup> in Eq. (31) and integrating over the volume of the conductor gives

$$\mathcal{U}\_{\text{surf}}(\mathbf{r};\mathbf{r}\_{0}) = \oint\_{\mathbf{r}=1}^{\infty} \sum\_{m=1}^{2n+1} \mathbf{D}\_{n}^{m}(\mathbf{r}\_{0}) (\Delta \boldsymbol{\varPsi}) \mathbf{E}\_{n}^{m}(\boldsymbol{\varvarmu}) \mathbf{E}\_{n}^{m}(\boldsymbol{\upnu}) d\boldsymbol{v}(\mathbf{\uphat{r}\_{0}}),\tag{40}$$

where Helmholtz's decomposition has been used. We compute the above integral analytically by expanding <sup>Ψ</sup>ðr0<sup>Þ</sup> in terms of the product of Lamè functions E<sup>m</sup> <sup>n</sup> <sup>ð</sup>μÞE<sup>m</sup> <sup>n</sup> ðνÞ, namely

$$\Psi(\mathbf{r}) = \sum\_{n=1}^{\infty} \sum\_{m=1}^{2n+1} \psi\_n^m(\rho) \mathbf{E}\_n^m(\mu) \mathbf{E}\_n^m(\nu). \tag{41}$$

After a series of cumbersome algebra, it can be shown that [31]

Mathematical Foundation of Electroencephalography http://dx.doi.org/10.5772/68021 49

$$\mathcal{U}I\_{\text{surf}} = \sum\_{n=1}^{\approx} \sum\_{m=1}^{2n+1} \mathbf{D}\_n^m(\mathbf{r}\_0) I\_n^m \to\_n^m(\mu) \mathbf{E}\_n^m(\nu), \tag{42}$$

where

As a result, we expand the coefficients B<sup>m</sup>

<sup>n</sup> as

gm

, y<sup>0</sup> <sup>¼</sup> <sup>1</sup> 2h1h<sup>3</sup>

<sup>n</sup> <sup>¼</sup> <sup>Q</sup> � <sup>∇</sup><sup>r</sup>0E<sup>m</sup>

which incorporates characteristics, such as the geometrical and physical properties of the conductor, as well as the EEG measurements, the algebraic manipulations leading to the solution are somewhat a little more painless. Solving a linear system of six equations with six unknowns, it can be shown [34] that the dipoles position r<sup>0</sup> and moment Q depends only on

> g3 2 g1 1 þ g5 2 g3 1 � g2 1g4 2 g1 1g3 1

Qm <sup>¼</sup> hm h1h2h<sup>3</sup>

<sup>2</sup>Þð<sup>Λ</sup> � <sup>a</sup><sup>2</sup>

<sup>2</sup>ÞðΛ<sup>0</sup> � <sup>a</sup><sup>2</sup>

where Λ, Λ<sup>0</sup> are ellipsoidal parameters depending on the semi-axes ai, i ¼ 1,2,3.

X∞ n¼1

by expanding <sup>Ψ</sup>ðr0<sup>Þ</sup> in terms of the product of Lamè functions E<sup>m</sup>

After a series of cumbersome algebra, it can be shown that [31]

<sup>Ψ</sup>ðrÞ ¼ <sup>X</sup><sup>∞</sup>

n¼1

2 <sup>X</sup><sup>n</sup>þ<sup>1</sup> m¼1

2 Xnþ1 m¼1

coordinates, namely for n ¼ 2 and m ¼ 1,2, we obtain two constrains similar to Eq. (25), which

<sup>3</sup><sup>Þ</sup> <sup>Q</sup>1x<sup>0</sup> <sup>Λ</sup> � <sup>a</sup><sup>2</sup> 1

<sup>3</sup><sup>Þ</sup> <sup>Q</sup>1x<sup>0</sup> <sup>Λ</sup><sup>0</sup> � <sup>a</sup><sup>2</sup> 1

<sup>n</sup> <sup>ð</sup>r0ÞðΔΨÞE<sup>m</sup>

In the case of a continuously distributed current, we follow the example laid out for the sphere. Replacing Q by the primary current J<sup>p</sup> in Eq. (31) and integrating over the volume of the

D<sup>m</sup>

where Helmholtz's decomposition has been used. We compute the above integral analytically

ψ<sup>m</sup> <sup>n</sup> <sup>ð</sup>ρÞE<sup>m</sup> <sup>þ</sup> <sup>Q</sup>2y<sup>0</sup> <sup>Λ</sup> � <sup>a</sup><sup>2</sup> 2

<sup>þ</sup> <sup>Q</sup>2y<sup>0</sup> <sup>Λ</sup><sup>0</sup> � <sup>a</sup><sup>2</sup> 2

<sup>n</sup> <sup>ð</sup>μÞE<sup>m</sup>

<sup>n</sup> <sup>ð</sup>μÞE<sup>m</sup>

<sup>n</sup> <sup>ð</sup>μÞE<sup>m</sup>

� �

� �

for every <sup>m</sup> <sup>¼</sup> 1,2,3, respectively. Expressing two of the triple products <sup>E</sup><sup>m</sup>

1, g<sup>2</sup> 1, g<sup>3</sup> 1, g<sup>3</sup> 2, g<sup>4</sup> 2, g<sup>5</sup> <sup>2</sup> as

<sup>1</sup>Þð<sup>Λ</sup> � <sup>a</sup><sup>2</sup>

<sup>1</sup>ÞðΛ<sup>0</sup> � <sup>a</sup><sup>2</sup>

� �

gm

<sup>n</sup> <sup>ð</sup>r0Þ ¼ <sup>d</sup> <sup>B</sup><sup>m</sup>

Defining a new parameter gm

48 Electroencephalography

g3 2 g2 1 þ g4 2 g3 1 � g1 1g5 2 g2 1g3 1

are satisfied by the constants g<sup>1</sup>

g1

g2

<sup>2</sup> <sup>¼</sup> <sup>2</sup>ð<sup>Λ</sup> � <sup>a</sup><sup>2</sup>

<sup>2</sup> <sup>¼</sup> <sup>2</sup>ðΛ<sup>0</sup> � <sup>a</sup><sup>2</sup>

<sup>U</sup>surfðr;r0Þ ¼ ∮

� �

g1 1, g<sup>2</sup> 1, g<sup>3</sup> 1, g<sup>3</sup> 2, g<sup>4</sup> 2, g<sup>5</sup> <sup>2</sup> as

and

and

conductor gives

<sup>x</sup><sup>0</sup> <sup>¼</sup> <sup>1</sup> 2h2h<sup>3</sup> <sup>n</sup> for n ¼ 1 and 2 yielding eight relations in total.

, z<sup>0</sup> <sup>¼</sup> <sup>1</sup> 2h1h<sup>2</sup>

<sup>n</sup> ð35Þ

g4 2 g1 1 þ g5 2 g2 1 � g3 1g3 2 g1 1g2 1

<sup>1</sup> , ð37Þ

<sup>þ</sup> <sup>Q</sup>3z<sup>0</sup> <sup>Λ</sup> � <sup>a</sup><sup>2</sup> 3

<sup>þ</sup> <sup>Q</sup>3z<sup>0</sup> <sup>Λ</sup><sup>0</sup> � <sup>a</sup><sup>2</sup> 3

� �

<sup>n</sup> ðrÞ in Cartesian

, ð39Þ

<sup>n</sup> ðνÞdvð^r0Þ, ð40Þ

<sup>n</sup> ðνÞ, namely

<sup>n</sup> ðνÞ: ð41Þ

ð36Þ

ð38Þ

$$\begin{split} I\_{n}^{m} &= \int\_{h\_{2}}^{a\_{1}} \psi\_{n}^{m}(\rho) \frac{\mathbf{d}}{\mathbf{d}\rho} \left( \sqrt{\rho^{2} - h\_{3}^{2}} \sqrt{\rho^{2} - h\_{2}^{2}} \frac{\mathbf{d}}{\mathbf{d}\rho} E\_{n}^{m}(\rho) \right) \mathbf{d}\rho \\\\ &+ \sqrt{a\_{1}^{2} - h\_{3}^{2}} \sqrt{a\_{1}^{2} - h\_{2}^{2}} \Big( E\_{n}^{m}(a\_{1}) \dot{\psi}\_{n}^{m} - \dot{E}\_{n}^{m}(a\_{1}) \psi\_{n}^{m} \Big). \end{split} \tag{43}$$

As for the spherical conductor, additional information is required in order to evaluate the coefficients ψ<sup>m</sup> <sup>n</sup> : A possible approach to uniquely determine the coefficients <sup>ψ</sup><sup>m</sup> <sup>n</sup> is to consider some form of minimizer. For example, a widely used concept in medical imaging is the minimum principle, the assumption that the current should have minimum 'strength', mathematically expressed as the minimization of the L<sup>2</sup> norm of Jp:

#### 4. Discussion

We presented a brief but concise introduction of the mathematical terminology associated with the modality of EEG. Having in mind medical and health professionals, we start from the very beginning, presenting step by step the physics and mathematical formulation behind EEG in a simple manner, keeping the mathematical notation to a minimum. The tools and techniques needed in order to derive at the presented results are intentionally not incorporated for two reasons. First of all, the procedure deriving at these formulas is not an easy task in general. Secondly, the main focus of this work is to display the beauty of the final expressions for every problem, showing how every single piece of information is encoded within these formulas and by what means the extraction of conclusions is accomplished.

Our introduction starts with the most elementary model possible, which, at the same time, is also the most straightforward and understandable of models. Representing the brain as a homogeneous sphere is an unrealistic assumption but serves an important task. At present, it is the only geometry for which the electromagnetic fields generated by a dipole source are exactly known in closed form. Further, it is needed in order to be able to draw conclusions when we move to build models of higher complexity. It therefore reveals disparities between forthcoming models, but more substantially, it allows us to test the reliability of the introduced algorithms in a straightforward and timely matter. Moreover, an analytic benchmark problem is provided, which can be used to test existing and new formulations.

For example, based on the homogeneous spherical model it is shown in Ref. [35] to what degree deformations present at the conductor's surface affect EEG measurements. Although the EEG data are evaluated in a view to a deformed conductor, the calculations are accomplished based on the spherical geometry, furnishing a fast analytic algorithm prone to almost minimum error. Another characteristic example of rigorous mathematical analysis is the quantitative description of the non-uniqueness for the EEG inverse problem, presented in Ref. [36]. Therein, splitting the current into components, the authors prove that none of those components contributes to both the electric and scalar magnetic potential; in other words, recordings of EEG and MEG do not contain overlapping information about the current. However, aforementioned property holds no longer true if the spherical conductor is disregarded.

Analogous conclusions are valid for the ellipsoidal geometry as well. For example, the authors consider in Ref. [37] the frequent case when clinical data of unknown origin are implemented in computational simulations. We mentioned earlier that there exist a plethora of combinations of the product a1, a2, a<sup>3</sup> furnishing the same value for the volume of the brain. If we look at the instance where EEG measurements originate from a brain with fixed values a1, a2, a<sup>3</sup> but are interpreted in the sequel by different values, what would be the error? Well, it turns out that the error can reach as high as 20%, depending on the position and strength of the primary current.

The error analysis presented in Ref. [37] can be considered as rather straightforward, since both ellipsoids under consideration were confocal, that is, members of the same ellipsoidal system enjoying the same foci. In plain words, no member of a confocal family touches another. Consequently, there exists a single curve, which cuts both ellipsoids normally, and the corresponding intersecting points consist of the most proximate pair between the two ellipsoids in each direction. These two points are employed in the analysis calculating the aforementioned error. But what would happen if the two ellipsoids would not be considered to be confocal? This interesting case is also more realistic.

In order to provide an answer to the latter, a sophisticated correspondence is needed connecting two points on the surface of two ellipsoids which are now non-confocal. This means that there probably exists a point shared by both ellipsoids. In Ref. [38], the authors investigated the effect implied by a deviation of the eccentricities of the ellipsoidal model on the electric potentials registered as the EEG data. In this case, the error reaches high values up to almost 100%.

Turning our attention from the geometrical deformation of the conductor model, to the physical assumption of homogeneity, we acknowledge the significance of the non-homogeneity imposed by the layers of different conductivity that cover the host tissue of the EEG source. The conductive elements that constitute the scalp, the scull and the meninges, which interfere between the EEG measurements and the cerebrum, are affected by the electromagnetic field produced by activation of the source. Hence, they induce a volume current that perturbs the total electric potential registered on the EEG receptors on the scalp. The effect of this physical perturbation of the potential has been studied by incorporating a layered conductivity profile in all the models discussed in the present review, by characterizing each layer by a distinct but constant conductivity value.

Indicatively, we infer that switching to a layered ellipsoidal model of the head-brain system, the functional form of the electric potential, is basically unaltered. One of the authors has showed [39] that the conductivity profile of the layered structure enters the potential formula by normalizing each term by a constant which incorporates the conductivity jumps across the interfaces and the geometrical characteristics of the layers. Analogous results show a similar effect on the inhomogeneous conductor [31]. Hence, the formula of the electric potential that will serve as the stepping stone for the inverse calculations is the one that corresponds to the most realistic inhomogeneous models that acknowledge the layered conductivity profile of the head-brain system. A promising challenge for future investigations refers to incorporating the anisotropic conductivity profile, where the conductivity varies with the direction into each separate compartment, modelling, for example, the different conductivity of the white matter of the brain than that of grey matter.

Using tensor conductivity for modelling the brain anisotropy is one of the many analytical mathematical challenges of this fascinating area of functional brain imaging while creative mathematical modelling has a lot to contribute in the close examination of EEG theory and applications.
