Mathematical Foundation of Electroencephalography

Abstract

Electroencephalography (EEG) has evolved over the years to be one of the primary diagnostic technologies providing information concerning the dynamics of spontaneous and stimulated electrical brain activity. The core question of EEG is to acquire the precise location and strength of the sources inside the human brain by knowledge of an electrical potential measured on the scalp. But in what way is the source recovered? Leaving aside the biological mechanisms on the cellular level responsible for the recorded EEG signals, we pay attention to the mathematical aspects of the narrative. Our goal is to provide a brief and concise introduction of the mathematical terminology associated with the modality of EEG. We start from the very beginning, presenting step by step the mathematical formulation behind EEG in a simple and clear manner, keeping the mathematical notation to a minimum. Whilst we serve only the key relations for the described problems, we focus specifically on the limitations of each modelling approach. In this fashion, the reader can appreciate the beauty of the formulas presented and discover every single piece of information encoded within these formulas.

1. Introduction

The human brain is a remarkable and fascinating organ exhibiting a tremendous complexity. It makes us unique and defines who we are. In spite of our scientific and technological progression, we do not know the particulars of its operating, and as we delve into its secrets, various surprises emerge, for example, nearly 100 previously unidentified brain areas have been recently discovered [1]. Consisting of an inconceivable network of interconnected nerve cells and fibres, continuously transporting and processing information, the brain is extremely vulnerable and requires paramount protection. Several layers of safety are incorporated starting with three connective sheets of tissue, called the meninges containing the cerebrospinal fluid, followed by plates of bones, the skull. Above machinery safeguards the brain from mechanical damage. On the other hand, a sophisticated barrier within the brain provides a natural defence against toxic or infective agents.

Examining the anatomy of the brain, we recognize three distinct regions. The largest part of the brain is the cerebrum, divided into two hemispheres. The outermost layer of the cerebrum is the cerebral cortex, consisting of four lobes. Cognitive awareness has its origins here. The second largest part of the brain is the cerebellum, located underneath the cerebrum and responsible for motor control and learning. Last, but not least, an integral part of the brain connecting the cerebrum with the spinal cord is the brainstem, regulating reflexes and crucial, basic life functions. Detailed information can be found in [2].

The operational status of the brain is based on an alternating chain of electrical and chemical events. On the microscopic level, encoding and transmitting of information via electrochemical signals is achieved by the active participation of neuronal and non‐neuronal constituents. Brain cells communicate through synaptic transmissions by controlling chemical transmitters or ionic currents which flow across their membranes. As a consequence, an electromagnetic field is generated. For a far‐reaching introduction on the subject, see [3]. The question at hand now lies in the possibility exploiting these provoked fields. It seems only reasonable that if a substantial number of cells form a critical mass, which activates synchronously, the emerging electric and magnetic field should be detectable. This is indeed the case and a deeply rooted concept in electrophysiology [4, 5]. From an electrofunctional point of view, the ionic micro‐flow within a single brain cell creates an opposite polarity between two point electrical charges very close together, leading to the notion of a dipole, a physical quantity one could say consists of the ‘fundamental unit’, which produces the observed fields. Dipoles are characterized by a vector called moment, the product of the charge and distance, visualized as an arrow pointing from a minimum (negative charge) to a maximum (positive charge), ergo featuring direction and magnitude. For that reason, it may be argued that the macroscopic description of the brain's activity is best achieved when simulated as an array of dipoles, that is, a non‐uniform distribution of positive and negative charges. According to the latter, if a small neighbourhood is stimulated, an excessive number of dipoles concur and their electric fields would add or cancel one another depending on the direction. This complicated and difficult situation can be avoided by introducing the concept of the equivalent current dipole (ECD), namely a single dipole which generates the identical electric field as all of the individual dipoles together, hence summarizing the net effect of all microscopic currents located in the distinct region of the brain under consideration. This is a widely used approximation concept in the framework of neuroelectromagnetism [6, 7]. On the other hand, when the exertion is no longer confined to a focal region of the brain, then every one of these regions is simulated by an equivalent current dipole, leading to a distribution of sources.

The main task and problem is to correlate active regions with associated generated electric fields. This essential step is closely connected with the installation of physical structures, namely a boundary or number of boundaries enclosing distinct regions with specific physical characteristics, such as conductivity. The head model obtained is termed the volume conductor model. Clearly, the level of details incorporated into the head model provides an analogous degree of operational freedom when it comes to investigate how the fields generated by brain cells are transmitted through various biological tissues towards measurement apparatus. As a result, the volume conductor model consists of the physical foundation for source analysis, which is categorized into two major problems. The first one is associated with the calculation of the electric potential, generated by known electrochemical sources within the brain, at precise points at the surface of the scalp. This is the forward electroencephalographic (EEG) problem [8–11]. The forward EEG problem has been extensively scrutinized for over 60 years since Wilson and Bayley [12] attempted to quantify the interplay between neuronal activity and the potentials they generate at the scalp. The reconstruction of the sources responsible for the recorded values is called the corresponding inverse EEG problem [13].

As of today, a high level of details can only be achieved with the aid of numerical models, which are generally categorized into boundary element models (BEMs) and finite element models (FEMs). Whereas boundary element models are adequate to portray major tissue compartments, such as the cerebrum and skull, they fail to represent detailed anatomical information within the compartments, such as the cerebral folding [14, 15]. Finite element models, on the other hand, are efficient in capturing these details, but are labour intensive and computationally demanding [16, 17].

Nonetheless, in order to gain a deeper comprehension of the problem a rigorous mathematical analysis is essential in providing a vital step towards the recognition of the underlying phenomena as well as identifying the limitations of the developed algorithms. The importance of mathematical analysis cannot be emphasized enough, since (i) it allows testing the impact of modifications regarding various variables upon the output of the system and provides further insight into underlying physical behaviour. (ii) It serves as validation tools for the numerical models.

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 E+\frac{\partial B}{\partial t}=0,\nabla \times H=J+\frac{\partial D}{\partial t},\nabla \cdot B=0,\nabla \cdot D=\rho .$E1

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 ($\cdot$) and cross ($\times$) 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 $\rho$ 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

$D=\epsilon \mathbf{\text{E}},H={\mu }^{-1}B,$E2

where $\epsilon$ is the relative permittivity (dielectric constant) of the material and shows how strong the material influences the electric field $E.$ Similarly, $\mu$ 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 E=0,\nabla \times H=J,\nabla \cdot B=0,\nabla \cdot D=\rho .$E3

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

$\nabla \times E=0,\nabla \times B=\mu J,\nabla \cdot B=0,\nabla \cdot D={\epsilon }^{-1}\rho .$E4

Let us focus on the first of Eq. (4). The curl of a vector, that is, the operator $\nabla \times$ 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,

$E=-\nabla U.$E5

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 J=0$E6

and the current density of $J$ is said to be solenoidal and expresses the steady‐state condition that the charge density $\rho$ is not changing in time. Moreover, an applied field in a resistive material, such as the brain, will induce a current of density $J_{\text{i}}$, 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]),

$J_{\text{i}}=\sigma E,$E7

where the proportionality constant $\sigma$ is termed conductivity. Therefore, the total current is given as the sum of the primary current $J_{\text{p}}$ responsible for the electric field, and the induced, or secondary, current $J_{\text{i}}$ given by Eq. (8), as

$J=J_{\text{p}}+\sigma E.$E8

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

$\nabla \cdot \left(\sigma \nabla U\right)=\nabla \cdot J_{\text{p}},$E9

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

$\nabla \cdot \left(\sigma \nabla U\right)=\nabla U\cdot \nabla \sigma +\sigma \Delta U,$E10

where the symbol $\Delta$ is called the Laplacian operator.

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 $\nabla \sigma$ 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_{\text{p}}$ by Poisson's equation,

$\Delta U={\sigma }^{-1}\nabla \cdot J_{\text{p}}.$E11

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_{\text{p}}$ multiplied by 1 over $\sigma$. In view of Eq. (11), the forward EEG problem is now formulated as follows. Given the primary current density $J_{\text{p}}$, calculate the electric potential on the surface of the medium. On the contrary, calculating $J_{\text{p}}$ 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 $V_1,V_2$ the wave and its normal derivative must be continuous at the interface, namely

$U_1=U_2\text{and}{\sigma }_1{\partial }_{\nu }{U}_1={\sigma }_2{\partial }_{\nu }{U}_2\text{on}S,$E12

where ${\sigma }_1,{\sigma }_2$ denote the conductivity of $V_1,V_2$, 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 ${\partial }_{\nu }U$ on $S$ must cancel out, that is,

${\displaystyle \oint {\partial }_{\nu }U\text{ d}S=0}.$E13

Further note that, by virtue of ${\partial }_{\nu }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_0$ and moment $Q$, we consider

$J_{\text{p}}=Q\delta \left(r-r_{\text{0}}\right),$E14

where the functional $\delta \left(r\right)$ is the Dirac measure, a concept included in the analysis by the reason of representing the concentration of $J_{\text{p}}$ to a single point [21].

3.1. The homogeneous spherical brain

The simplest possible geometry in order to represent the brain consists of a spherical homogeneous conductor with radius $a$ and conductivity $\sigma$. 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 U={\sigma }^{-1}Q\cdot \nabla \delta \left(r-r_0\right),r<a$E15

supplemented by the condition

${\partial }_{r}U=0,r=a,$E16

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.¹ Note that the compatibility condition (13) is automatically satisfied by Eq. (16).

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

$U_{\text{surf}}\left(\widehat{\mathbf{r}},r_0\right)={\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=-n}^n{\text{A}}_n^m\left({\widehat{\mathbf{r}}}_0\right)Y_n^m\left(\widehat{\mathbf{r}}\right)}},{\text{ A}}_n^m\left({\widehat{\mathbf{r}}}_0\right)=\frac{Q}{\sigma n{a}^{n+1}}\cdot {\nabla }_{r_0}\left(r_0^n{\overline{Y}}_n^m\left({\widehat{\mathbf{r}}}_0\right)\right),$E17

where a hat ‘^’ denotes the unit vector, $Y$ and $\overline{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:

$U_{\text{surf}}\left(\widehat{\mathbf{r}},r_0\right)=\frac{Q}{4\pi \sigma }\cdot \left(2\frac{R}{R^3}+\frac{1}{aR}\frac{R\widehat{\mathbf{r}}+R}{R+\widehat{\mathbf{r}}\cdot \mathbf{R}}\right),R=a\widehat{\mathbf{r}}-r_0,$E18

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 $U_{\text{surf}},$ given via Eq. (17), from which we have to calculate the position $r_{\text{0}}=\left(x_0,y_0,z_0\right)$ and moment $Q=\left(Q_x,Q_y,Q_z\right)$ of the dipole. In our case, this is easily achieved, noting that every detail regarding $r_0$ and $Q$ are encoded into the coefficients ${\text{A}}_n^m$ of expansion (17). The latter are evaluated with the aid of the orthogonality condition

${\text{A}}_n^m={\displaystyle \oint U_{\text{surf}}\left(\widehat{\mathbf{r}},r_0\right){\overline{Y}}_n^m\left(\widehat{\mathbf{r}}\right)\text{d}S\left(\widehat{\mathbf{r}}\right)}.$E19

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 $U_{\text{surf}}$ 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 ${\text{A}}_n^m$, given via Eq. (17), for n = 1 and 2 providing eight relations in total. The first three of them, namely ${\text{A}}_1^{-1},{\text{ A}}_1^0$ and ${\text{A}}_1^1$, are proportional to the components $Q_x,Q_y,Q_z$ of $Q,$ and thus we find

$Q_x=\sqrt{\frac{2\pi }{3}}\sigma a^2\left({\text{A}}_1^{-1}-{\text{A}}_1^1\right),Q_y=-i\sqrt{\frac{2\pi }{3}}\sigma a^2\left({\text{A}}_1^{-1}+{\text{A}}_1^1\right),Q_z=2\sqrt{\frac{\pi }{3}}\sigma a^2{\text{A}}_1^0,$E20

where $i$ denotes the imaginary unit. It is not hard to show that whereas the difference ${\text{A}}_1^{-1}-{\text{A}}_1^1$ is real valued, the corresponding sum is imaginary. The remaining five coefficients ${\text{A}}_2^m$ 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{{\text{A}}_2^{-2}}{{\text{A}}_1^{-1}}-\frac{{\text{A}}_2^2}{{\text{A}}_1^1}\right),y_0=-i\frac{a}{\sqrt{5}}\left(\frac{{\text{A}}_2^{-2}}{{\text{A}}_1^{-1}}+\frac{{\text{A}}_2^2}{{\text{A}}_1^1}\right),z_0=\frac{2a}{\sqrt{5}}\frac{1}{{\text{A}}_1^1}\left({\text{A}}_2^1-\frac{{\text{A}}_2^2{\text{A}}_1^0}{\sqrt{2}{\text{A}}_1^1}\right).$E21

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]

$U_{\text{surf}}=A_1{x}_0+A_2{y}_0+A_3{z}_0+B_1{x}_0{}^2+B_2{y}_0{}^2+B_3{z}_0{}^2+C_1{x}_0{y}_0+C_2{x}_0{z}_0+C_3{y}_0{z}_0+\dots$E22

where again the first three coefficients $A_1,A_2,A_3$ 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 $Q_x,Q_y,Q_z$ into $B_1,B_2,B_3,C_1,C_2,C_3$, we arrive at the linear but overdetermined system $\mathbb{A}{r}_0=b$ with

$\mathbb{A}=\left[\begin{array}{ccc}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{array}\right],b=\frac{6a^2}{5}\left[\begin{array}{c}{B}_1\\ B_2\\ B_3\\ C_1\\ C_2\\ C_3\end{array}\right].$E23

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

$x_0=\frac{3}{2A_1{A}_2{A}_3}\left(A_{1}a\right)\cdot \tilde{\mathbf{b}},y_0=\frac{3}{2A_1{A}_2{A}_3}\left(A_{2}a\right)\cdot \left({ℝ}_z\tilde{\mathbf{b}}\right),y_0=\frac{3}{2A_1{A}_2{A}_3}\left(A_{3}a\right)\cdot \left({ℝ}_y\tilde{\mathbf{b}}\right),$E24

where $a={\left(A_1,A_2,A_3\right)}^T,\tilde{\mathbf{b}}={\left(-C_1,C_2,C_3\right)}^T,$$T$ denotes transposition, whereas ${ℝ}_y$ and ${ℝ}_z$ 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{array}{l}2A_1{A}_2{A}_3\left(B_3-B_2\right)+3A_1\left(A_2^2-A_3^2\right)C_2+3\left(A_2^2+A_3^2\right)\left(A_3{C}_1-A_2{C}_3\right)=0,\\ -2A_1{A}_2{A}_3\left(B_2+2B_3\right)+3A_2\left(A_3^2-A_1^2\right)C_3+3\left(A_1^2+A_3^2\right)\left(A_1{C}_2-A_3{C}_1\right)=0.\end{array}$E25

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_0={\mathbb{A}}^{+}b+\left(\mathbb{I}-{\mathbb{A}}^{+}\mathbb{A}\right)Y,$ where ${\mathbb{A}}^{+}$ is the pseudoinverse of $\mathbb{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 ${\mathbb{A}}^{+}\mathbb{A}=\mathbb{I}$ holds true, leading to the minimum norm solution $r_0={\mathbb{A}}^{+}b.$

However, in ‘reality’ EEG recordings provide values for the coefficients $A_n^m$, 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 $U_j$ corresponding to dipoles located at $r_{0j}$ and strength $Q_j.$ 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

$U_{\text{surf}}\left(\widehat{\mathbf{r}},r_0\right)=\frac{1}{4\pi \sigma }{\displaystyle \sum _{n=1}^{\infty }\frac{2n+1}{n{a}^{n+1}}\left(Q\cdot {\nabla }_{r_0}\right)r_0^n{P}_n\left(\widehat{\mathbf{r}}\cdot {\widehat{\mathbf{r}}}_0\right)},$E26

employing the addition theorem for Legendre functions $P_n.$ The same expression is valid on the occasion of N dipoles, that is,

$U_{\text{surf}}\left(\widehat{\mathbf{r}},r_{0j}\right)=\frac{1}{4\pi \sigma }{\displaystyle \sum _{n=1}^{\infty }\frac{2n+1}{n{a}^{n+1}}{\displaystyle \sum _{j=1}^N\left(Q_j\cdot {\nabla }_{r_{0j}}\right)r_{0j}^n{P}_n\left(\widehat{\mathbf{r}}\cdot {\widehat{\mathbf{r}}}_{0j}\right)}}.$E27

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

$\left(Q\cdot {\nabla }_{r_0}\right)r_0^n{P}_n\left(\widehat{\mathbf{r}}\cdot {\widehat{\mathbf{r}}}_0\right)={\displaystyle \sum _{j=1}^N\left(Q_j\cdot {\nabla }_{r_{0j}}\right)r_{0j}^n{P}_n\left(\widehat{\mathbf{r}}\cdot {\widehat{\mathbf{r}}}_{0j}\right)}.$E28

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

If we abandon the assumption that the primary current $J_{\text{p}}$ 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_0$ and oriented along an arbitrary direction $\widehat{\mathbf{\alpha }}.$ The primary current is then approximated as $J_{\text{p}}\simeq Q+tA,$ where t is a variable taking values in the interval [‐L,L], whereas $A$ is the directional derivative of $J_{\text{p}}$ along $\widehat{\mathbf{\alpha }}.$ Replacing the approximation of $J_{\text{p}}$ 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_0$ and perpendicularly oriented to the vector $r_0,$ the primary current can be approximated as $J_{\text{p}}\simeq Q+r\cdot \tilde{\mathbf{D}}.$ The quantity $\tilde{\mathbf{D}}$ is called dyadic, a second‐order tensor containing physical or geometrical information. The surface potential is provided replacing $J_{\text{p}}$ 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_{\text{p}}$ 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, $\nabla \Psi (r_0)$ and $A(r_0),$ respectively. As a result, assuming $J_{\text{p}}=\nabla \Psi +\nabla \times A$ under the condition that $\nabla \cdot A=0,$ we find

$U_{\text{surf}}\left(\widehat{\mathbf{r}},r_0\right)={\displaystyle \oint {\displaystyle \sum _{n=1}^{\infty }{C}_n^m\left(\Delta \Psi \right)r_0^n{P}_n\left(\widehat{\mathbf{r}}\cdot {\widehat{\mathbf{r}}}_0\right)\text{d}v({\widehat{\mathbf{r}}}_0)}}.$E29

Note that the electric potential depends only on $\Delta \Psi (r_0)$. The integral in Eq. (29), obtained using Gauss theorem and integration by parts, can be computed analytically if we expand $\Psi (r_0)$ in terms of spherical harmonics, namely

$\Psi \left(r_0\right)={\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=-n}^n{\psi }_n^m\left(r_0\right)Y_n^m\left({\widehat{\mathbf{r}}}_0\right)}},$E30

furnishing

$U_{\text{surf}}={\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=-n}^n\left({\alpha }_n\underset{n}{\overset{m}{\dot{\psi }}}-{\beta }_n{\psi }_n^m\right)Y_n^m\left(\widehat{\mathbf{r}}\right)}}.$E31

It is remarkable that the above coefficients can be determined only under the assumption that $\Psi (r_0)$ is a bi‐harmonic function, namely a solution to the bi‐harmonic operator $\Delta \Delta .$ However, the precise description of the coefficients ${\psi }_n^m$ 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 $a_1,a_2,a_3$, 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

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 $\rho$. This is secured if we use the boundary of our domain to be the so‐called reference ellipsoid and construct the ellipsoidal system that is based on it. Further, denote by $h_1,h_2,h_3$ the three semifocal distances from which the orthogonal ellipsoidal coordinate system $\text{(}\rho ,\nu ,\mu \text{)}$ is derived. Each confocal ellipsoidal surface is defined by a constant value of the ‘radial’ variable $\rho \in \text{(}{h}_2,\infty \text{)}$, with $\rho =a_1$ indicating the reference ellipsoid.

Consider in what follows a homogeneous ellipsoidal conductor with semi‐axes $a_i,$ i = 1,2,3 and conductivity $\sigma$ 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 U={\sigma }^{-1}Q\cdot \nabla \delta \left(r-r_0\right),h_2<\rho <a_1,$E32

${\partial }_{\rho }U=0,\rho =a_1.$E33

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

$U_{\text{surf}}\left(r;r_0\right)=\frac{1}{\sigma }{\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=1}^{2n+1}{\text{B}}_n^m\left(r_0\right){\text{E}}_n^m\left(\mu \right){\text{E}}_n^m\left(\nu \right)}},{\text{ B}}_n^m\left(r_0\right)=\frac{Q\cdot {\nabla }_{r_0}{\mathbb{E}}_n^m\left(r_0\right)}{d}.$E34

where ${\text{E}}_n^m$ are the Lamé functions, ${\mathbb{E}}_n^m\left(r\right)$ 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_0$ 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_0$ of the dipole is encoded into the coefficients of expansion (31). As a result, we expand the coefficients ${\text{B}}_n^m$ for n = 1 and 2 yielding eight relations in total. Defining a new parameter $g_n^m$ as

$g_n^m=Q\cdot {\nabla }_{r_0}{\mathbb{E}}_n^m\left(r_0\right)=d{\text{ B}}_n^m$E35

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_0$ and moment $Q$ depends only on $g_1^1,g_1^2,g_1^3,g_2^3,g_2^4,g_2^5$ as

$x_0=\frac{1}{2h_2{h}_3}\left(\frac{g_2^3}{g_1^2}+\frac{g_2^4}{g_1^3}-\frac{g_1^1{g}_2^5}{g_1^2{g}_1^3}\right),y_0=\frac{1}{2h_1{h}_3}\left(\frac{g_2^3}{g_1^1}+\frac{g_2^5}{g_1^3}-\frac{g_1^2{g}_2^4}{g_1^1{g}_1^3}\right),z_0=\frac{1}{2h_1{h}_2}\left(\frac{g_2^4}{g_1^1}+\frac{g_2^5}{g_1^2}-\frac{g_1^3{g}_2^3}{g_1^1{g}_1^2}\right)$E36

and

$Q_m=\frac{h_m}{h_1{h}_2{h}_3}{g}_1^m,$E37

for every m = 1,2,3, respectively. Expressing two of the triple products ${\mathbb{E}}_n^m\left(r\right)$ 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_1^1,g_1^2,g_1^3,g_2^3,g_2^4,g_2^5$ as

$g_2^1=2\left(\Lambda -a_1^2\right)\left(\Lambda -a_2^2\right)\left(\Lambda -a_3^2\right)\left(\frac{Q_1{x}_0}{\Lambda -a_1^2}+\frac{Q_2{y}_0}{\Lambda -a_2^2}+\frac{Q_3{z}_0}{\Lambda -a_3^2}\right)$E38

and

$g_2^2=2\left(\Lambda \text{'}-a_1^2\right)\left(\Lambda \text{'}-a_2^2\right)\left(\Lambda \text{'}-a_3^2\right)\left(\frac{Q_1{x}_0}{\Lambda \text{'}-a_1^2}+\frac{Q_2{y}_0}{\Lambda \text{'}-a_2^2}+\frac{Q_3{z}_0}{\Lambda \text{'}-a_3^2}\right),$E39

where $\Lambda ,{\Lambda }^{\prime }$ are ellipsoidal parameters depending on the semi‐axes $a_i,$ 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_{\text{p}}$ in Eq. (31) and integrating over the volume of the conductor gives

$U_{\text{surf}}\left(r;r_0\right)={\displaystyle \oint {\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=1}^{2n+1}{\text{D}}_n^m\left(r_0\right)\left(\Delta \Psi \right){\text{E}}_n^m\left(\mu \right){\text{E}}_n^m\left(\nu \right)dv\left({\widehat{\mathbf{r}}}_0\right)}}},$E40

where Helmholtz's decomposition has been used. We compute the above integral analytically by expanding $\Psi (r_0)$ in terms of the product of Lamè functions ${\text{E}}_n^m\left(\mu \right){\text{E}}_n^m\left(\nu \right),$ namely

$\Psi \text{(}r\text{)}={\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=1}^{2n+1}{\psi }_n^m\left(\rho \right){\text{E}}_n^m\left(\mu \right){\text{E}}_n^m\left(\nu \right)}}.$E41

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

$U_{\text{surf}}={\displaystyle \sum _{n=1}^{\infty }{\displaystyle \sum _{m=1}^{2n+1}{\text{D}}_n^m\left(r_0\right){{\rm I}}_n^m{\text{ E}}_n^m\left(\mu \right){\text{E}}_n^m\left(\nu \right)}},$E42

where

${{\rm I}}_n^m={\displaystyle \underset{h_2}{\overset{a_1}{\int }}{\psi }_n^m\left(\rho \right)\frac{\text{d}}{\text{d}\rho }\left(\sqrt{{\rho }^2-h_3^2}\sqrt{{\rho }^2-h_2^2}\frac{\text{d}}{\text{d}\rho }{E}_n^m\left(\rho \right)\right)\text{d}\rho }+\sqrt{a_1{}^2-h_3^2}\sqrt{a_1{}^2-h_2^2}\left(E_n^m\left(a_1\right){\dot{\psi }}_n^m-{\dot{E}}_n^m\left(a_1\right){\psi }_n^m\right).$E43

As for the spherical conductor, additional information is required in order to evaluate the coefficients ${\psi }_n^m.$ A possible approach to uniquely determine the coefficients ${\psi }_n^m$ 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_2$ norm of $J_{\text{p}}\text{.}$

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 $a_1,a_2,a_3$ 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 $a_1,a_2,a_3$ 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.