The model and the Inverse Electroencephalographic Problem

In the simple volume conductor model (studied in (^{Sarvas, 1987}) based on the results in (^{Geselowitz, 1967})), the brain is considered as conducting medium of electrical current, in which there is also a generating mechanism of other biological currents produced by neuronal activity, called impressed currents.

We will denote by Ω1 the region occupied by the brain, and by ∂Ω1=S1 its border, corresponding to the cortex. σ1 denotes the Ohmic current conductivity (in a normal brain is considered to be constant and equal to the one corresponding to the saturated salt water (^{Nunez & Srinivasan, 2006})) and μ denotes the magnetic permeability. From the perspective of a macro spatial scale, the EEG measurements detect average synaptic source activity (^{Nunez, Nunez, & Srinivasan, 2019}), so the constant conductivities reflect the average effect of microscopic spatial fluctuations. Additionally, we suppose that the electric field generated in the brain is due to their conductive properties as physical medium, and electrical sources originated by neuronal activity. Finally, u1 will be the electric potential in the brain volume Ω1. On the other hand, the rest of the head is considered also as a homogeneous conducting medium Ω2, with outer border ∂Ω2 =S2 corresponding to the scalp. Its average conductivity, considered constant, is σ2. In this region there are no sources of electrical activity (see Fig. 1). Analogously, u2 will be the electric potential distribution in Ω2. We set Ω=Ω1∪Ω2.

In this way, other layers such as the skull or the spinal brain fluid, among others, are not considered. The outer medium (outside the head) is supposed to have a vanishing conductivity. Hence, the simple quasi-static model describes the behavior of potentials u₁ and u₂, in the following way:

Where ∂ui∂nj denotes the normal derivative of u_i on S_j with respect to the normal unitary vector nj, outer to Ωj, i,j=1,2. Here, the forward problem consists on, given a source f, finding the potentials u₁ and u₂. Thus, the corresponding inverse problem consists on, given a measurement V, finding an equivalent source f satisfying ( 1 ) - ( 5 ) in such a way that the corresponding potential simultaneously fulfills the condition u2S2=V.

We consider the Hilbert spaces L2Ωi, L2Si and L2Ω of square summable functions defined on Ωi, Si and Ω; respectively i=1,2. Analogously, we will denote by H1Ωi and H1Ω the corresponding Sobolev spaces i=1,2. We will denote also by H1Si the subspaces of L2Si given by the traces in Si of functions in H1Ωi i=1,2.

The superindex (1) will indicate subspace corresponding to functions orthogonal to constants (with respect to the appropriate inner product). Namely, if W is a function Hilbert space with inner product ,W, then W1=w∈W:w,1W=0, where 1 is the indicator function of the corresponding domain. Also, the superindex ⊥ denotes the orthogonal complement in W of a set F in W, that is, F⊥=w∈W:w,fW=0, ∀f∈F.

In the neuronal activity source estimation problem, via model ( 1 ) - ( 5 ), the operative formulation requires to find an operator A which associates to each neuronal activity source f (in a certain class F) the measurement V. Namely, the inverse source estimation problem reduces to solve the operational equation

Where V=u2S2 is the potential distribution measurement on S₂.

Hence, solving the inverse source estimation problem means that, starting from the instantaneous electroencephalographic measurement V=u2S2, a source f “reproducing this measurement” could be find, i.e., substituting f in the model ( 1 ) - ( 5 ), there exists a unique solution u1f,u2f  satisfying u2fS2=V.

Given a class F of functions defined on the brain, the image of F by the operator A will be called admissible data set. This is exactly the set of measurements which could be reproduced by sources in the class F. Furthermore, the class F will be said to be a unicity class if operator A restricted to F is injective. In what follows the admissible data set associated to a sources class F will be denoted by MF:

In (^{Fraguela Collar, Oliveros Oliveros, Morín Castillo, & Conde Mones, 2015}) was proved that operator A is compact from L21Ω1 to L21S2. In addition, in Theorem 2.3 was shown that KerA⊥=H1Ω1, where HΩ1 is the closed subspace in L2Ω1 of harmonic functions. From this result yield the following conclusions, which will be important in what follows:

We note that another unicity class for the source estimation problem was obtained in (^{Fraguela Collar, Oliveros Oliveros, Morín Castillo, & Conde Mones, 2015}) consisting on certain kind of piecewise constant sources.

Operational formulation of the inverse problem

In order to solve the source estimation problem for certain unicity classes in the context of the volume conductor model ( 1 ) - ( 5 ), we need to reduce the inverse problem to an equivalent operational formulation. In general, solving this inverse problem in a realistic geometry is a quite involved mathematical problem. For the sake of simplicity, in order to explain our methodology more easily, we will consider a simple geometric model for the head.

We consider two concentric spheres. The interior of the inner sphere S₁, of radius R1>0 corresponds to the brain volume Ω1, and the interior of the outer sphere S₂, of radius R2>R1, corresponds to the whole head Ω and outlines the region Ω2 corresponding to the spherical crown (see Fig. 1 and Fig. 2).

This simple spherical model allows us to build explicit analytical expressions for the solutions of the volume conductor model ( 1 ) - ( 5 ) in terms of the Fourier series with respect to classical orthonormales bases, which eases the qualitative analysis, and also constitutes the basis of an algorithm for the numerical resolution of the forward and inverse problems.

We start from assuming that the source is in a unicity class F in L21Ω1. Note that, in L21Ω1, the forward problem ( 1 ) - ( 5 ) has an unique solution in the Sobolev space H11Ω, and hence its restriction to the border S₂ could be considered in the sense of Sobolev traces on S₂, (^{Mijailov, 1978}).

Next, we introduce some contour problems, which are auxiliary in solving the main inverse problem of source estimation, and associated to the corresponding solution. Note that these contour problems have no physical meaning; we use them for analyzing the operational formulation of the inverse problem. Their solutions will be interpreted in a weak sense. Fix a source f∈L21Ω1 and ψ∈L21S1.

Definition 1. The functions w1,v1∈H1Ω1, v2∈H1Ω2, satisfying

will be called weak solutions of problems ( 8 ), ( 9 ) and ( 10 ), respectively.

The following result holds.

Theorem 1. Problems ( 8 ), ( 9 ) and ( 10 ) have weak solution if and only if f∈L21Ω1 and ψ∈L21S1 . In this case, the weak solutions are unique in their respective spaces H11Ω1 and H11Ω2, and the following inequalities hold:

Where constants C₁, K₁ and C₂ do not depend on f and ψ, (see ^{Mijailov, 1978}).

Consequently, the following operators are well defined by using the weak solutions of problems ( 8 ), ( 9 ) and ( 10 ):

Making use of previous results and the compactness of the trace operators from H1Ω1 to L2S1, and from H1Ω2 to L2S1 and to L2S2, we conclude the compactness of the above operators, A₀, B₀, C₀ and D₀. It is convenient to consider the projection operators Pi:L2Si⟶L21Si, defined by:

Where Si denotes the Lebesgue measure of S_i. These averages could be considered as reference potentials, and these substractions are required in order to assure existence of the corresponding solutions. Thus we define the following operators:

By using these operators A,B,C and D, the solution of the inverse problem associated to the contour problem ( 1 ) - ( 5 ) can be obtained by solving the following system of operational equations:

Actually, equation ( 27 ) is equivalent to the Cauchy problem in region Ω2 for the Laplace operator with respect to the Cauchy data ∂u2∂n2=0 on S₂ and u2S2=V , where V∈L21S2. Once ψ is known, it can be substituted in ( 26 ), obtaining the operational equation:

Thus, equation ( 28 ) give us f, so system ( 26 ) - ( 27 ) is equivalent to inverse source estimation problem. Finally, the operator which relates the source f with the measurement V is given by:

The following result can be found in (^{Fraguela Collar, Oliveros Oliveros, Morín Castillo, & Conde Mones, 2015}).

Theorem 2. One has

This theorem justifies remarks a) and b) above. An important corollary follows also from it: if f∈F reproduces a measurement V, and the harmonic function h₀ is the unique harmonic source on Ω1 which also reproduces V, then h₀ is the orthogonal projection of f on H1Ω1, independently of the chosen class F, and f-h0 is the component of f in Ker A.

Hence, for the class F to be an unicity class is necessary and sufficient that

Where F-F=f1-f2 : f1,f2∈F and 0 is the vanishing function. Note that the class F need not to be a linear space. Indeed, F is a unicity class if A (or A) restricted to F is injective. This happens if for any f1,f2∈F with Af1=Af2 one gets f1=f2, but Af1-f2=Af1-Af2, so Af1=Af2 is equivalent to f1-f2∈KerA, which has to be trivial.

Admisible data method

The “Admissible Data Method” (ADM), as a regularization strategy (^{Hernandez Montero, Fraguela Collar, & Henry, 2019}), can be applied to any sources unicity class F. Concretely, in the context of the inverse source estimation problem linked to the instantaneous contour problem ( 1 ) - ( 5 ), when the admissible data set MF is contained into the admissible data set corresponding to the harmonic sources H1Ω1, consists of the following steps:

This methodology will be illustrated in the last section, in the case of the sources supported and harmonic in a neighborhood of the cortex.

Reduction of the instantaneous source estimation problem to a Cauchy data problem on cortex

In order to find the restriction to S₁ of the solution u₂ on Ω2, we will study the following contour problem, which is a Cauchy problem for the Laplace equation on Ω2, with vanishing Neumann contour condition and Dirichlet data V:

Where V is the potential distribution electroencephalographic measurement V in S₂ (the whole scalp) and u₂ is the electric potential in Ω2, at a given time instant. We proceed by computing formally the solution and studying the convergence of the obtained series (^{Mijailov, 1978}).

We consider the normalized spherical harmonics (^{Tijonov & Samarsky, 1980, pp. 765-778}):

Where Pnmθ,ϕ are the associated Legendre polynomials, and θ∈0,π, ϕ∈0,2π,1≤n<∞,-n≤ m≤n. Thus, V is given by

Where V_nm are the Fourier coefficients of V. It can be checked that the corresponding solution to problem ( 30 ) - ( 32 ) is given by

With convergence in L2Ω2. The solution of the Cauchy problem corresponds to restrict u₂ to S₁, which means to evaluate ( 35 ) at r=R₁, obtaining

In order to assure that the above expression for u2R1,θ,ϕ to be the trace in S₁ of a function in H1Ω1, the series composed with the squares of Fourier coefficients of ( 36 ), multiplied by n, must converge (^{Mijailov, 1978, p. 221}). Then, by using an elementary result on convergence of numerical series¹ we obtain the following condition on the Fourier coefficients of V:

Next we find out what conditions on V must be imposed in order for ψ=σ2σ1∂u2∂rr=R1 to be L21S1. ² It is easy to see that

Consequently, ( 38 ) must converge in L₂(S₁). Proceeding as before, we get the following result:

Theorem 3. One has ψ∈L21S1 if and only if the Fourier coefficients V_nm satisfy

All of this machinery allows us to convert the source estimation problem in the head into a source estimation problem in the brain from Cauchy data on the cortex. From now on we will suppose that the Fourier coefficients V_nm of V satisfy condition ( 39 ). From an operational point of view, it makes sense to solve equation ( 28 ) from ψ given by ( 38 ). That is

Which could be viewed as the identification of the source f, which reproduces the potential g on the cortex S₁, as it would be artificially measured.

In order to obtain an explicit expression for g=Cψ-Bψ, we make use of ( 38 ). Since operators B and C linear and continuous, it is enough to compute B(Ynmθ,ϕ) and C(Ynmθ,ϕ). For this purpose, in view of ( 18 ) and ( 19 ), we solve problems ( 9 ) and ( 10 ) with contour data Ynmθ,ϕ.

For B(Ynmθ,ϕ), we look for the solution v1 of the contour problem ( 9 ) of the form

And for C(Ynmθ,ϕ), we seek the solution v2 of the contour problem ( 10 ) of the form

In this way we get

From this expressions one easily gets

In summary, we conclude that, if the Fourier coefficients V_nm of measurement V satisfy condition ( 39 ), then the source estimation problem turns to solve equation Af = g, where g is given by ( 45 ). This problem can be interpreted as a source estimation problem in the brain, assuming it is isolated, and starting from the “measurement” g on the cortex. In what follows, we will apply this conclusion.

Identification of harmonic sources in the brain

In this section we consider a given instantaneous potential distribution V on the scalp, and we study under what conditions there exists an (unique) harmonic source f in the brain which reproduces V. By the way, we will find an explicit expression for f.

We start from the fact that any harmonic source f in Ω1 should take the following form:

Which is a development with respect to the orthonormal basis of spherical harmonics 2n+3R12n+3rnYnmθ,ϕ. Hence, the convergence of this series (in the L2Ω1 sense) is equivalent to the condition

Since operator A is linear and continuous, we get

In order to compute

In view of ( 17 ), it is required to solve the associated problem

Thus we look for a solution of the form

Where vr satisfies the equation

The general solution of this equation takes the form

Where v0 is solution for the associated homogeneous equation and vp is a particular solution. We propose a particular solution of the form

From ( 53 ) we obtain

Thus we get the condition k=n+2, so

And

The function v0 satisfies the homogeneous equation

And must be bounded at 𝑟=0, so takes the form

From the contour condition ( 51 ) we get

So

In this way we have finally found that

From ( 40 ) we must equate expressions ( 63 ) and ( 45 ). From the unicity of the coefficients, we obtain the Fourier coefficients F_nm in terms of the corresponding V_nm. Concretely, the final expression for f is:

Where the convergence must be in L2Ω1. Comparing expressions ( 48 ) and ( 64 ) we conclude the following result.

Theorem 4. There exists a biunivocal correspondence between the set of harmonic sources f in the brain Ω1 and the set of measurements V on the scalp (in L2S2) whose Fourier coefficients V_nm with respect to the orthonormalized spherical harmonics Ynmθ,ϕ satisfy:

Furthermore, given a measurement V satisfying condition ( 65 ), the unique harmonic source in Ω1 which reproduces it is given by ( 64 ).

We recall that from the previous theorem it follows that H1Ω1 is a sources unicity class, and the admissible data set corresponding to this class of harmonic sources orthogonal to the constants, denoted by MH1Ω1 following notation ( 7 ), turns to be the set of measurements V satisfying ( 65 ).

Note that, from Theorem 2 and 4 it follows that for any sources class F in L2Ω1 the corresponding admissible data set MF is contained into the set of measurements V in L2S2 which satisfy condition ( 65 ). Note also that this condition ( 65 ) is a quite strong smooth requirement on the spatial distribution of the measurements V, for each time instant.

Importance of the class of harmonic sources in the sources identification methodology on arbitrary unicity classes. The particular case of harmonic sources in a neighbourhood of the cortex

In this section we assume as before a given instantaneous potential distribution V on the scalp, and we study the existence of a source in the brain reproducing it. For this source f we assume it is supported and harmonic in a neigbourhood of the cortex. We stablish conditions on V for its existence, and we will find the admissible data set for this class.

Let us first clarify what will be understood for “a neighbourhood of the cortex” (see Fig. 3): Fix R ₀ satisfying 0<R0<R1. Ω10=0<r≤R0 represents the subcortical part of the brain, Ω11=R0<r≤R1 stands for a certain spherical neighbourhood of the cortex, and Ω2=R1<r≤R2 is the rest of the head as before. Analogously, the sphere S0=r=R0 is the border of Ω10, the sphere S1={r=R1} stands for the cortex and is the common border of Ω10 and Ω11, and S2={r=R2} stands for the scalp, where the measurements are taken, which is the outer border of Ω2.

It will be convenient to denote

We consider the class H1Ω11 of sources in the brain which are supported and harmonic in Ω1 (vanishing outside). It is easy to see that this sources f can be represented in the form (^{Tijonov & Samarsky, 1980, p. 772}):

Where the corresponding series of Fourier coefficients

should converge.

Next, we solve the auxiliary contour problem ( 8 ), which in view of ( 17 ) means compute A( f ). Since operator A is linear and continuous, it is enough to compute Ar^nYnmθ,ϕ and Ar^-n+1Ynmθ,ϕ. In order to find:

it is necessary to solve the associated contour problem

whose solution is imposed to take the form

It is easy to see that a (r) satisfies the equation

If we denote by a0r, w10 and a1r, w11 the corresponding restrictions of a (r) and w ₁ to Ω10 and Ω11, respectively, and consider the compatibility conditions on S ₀,

we arrive to functions a0r and a1r satisfy

and the contour conditions at r=R0

Also, the following additional conditions must be added:

Solving this problem for a0r and a1r gives us

Thus, we finally get

In a similar fashion one obtains

Equipped with ( 84 ) and ( 85 ), we can easily write the resulting formula for A ( f ):

Finally, we must equate (86) and ( 45 ) and apply the coefficients unicity in order to obtain

Let us define αn, βn, and γn by

It is easy to see that

Let us write ( 87 ) in the form

or equivalently

This yields

Finally, have the sum for any positive integer n and m from −n to n. If we take into account that f is harmonic in Ω1, we have the following conclusions:

In this way, we arrive to the equivalent condition:

which turns to be the necessary and sufficient condition previously obtained for V to be reproducible by an harmonic source defined in Ω1. That is, if there is an harmonic source f∈H1Ω11 which reproduces V, then ( 97 ) is fulfilled. In other words, ( 97 ) is a necessary condition.

Next, we will characterize the admissible data set corresponding to harmonic sources in a neighborhood of the cortex. By Theorem 2 we have that, if f∈H1Ω11 and h0∈H1Ω1 both reproduce V, then f-h0∈H1Ω1⊥, so

for n≥1, -n≤m≤n, where , denotes the inner product in L2Ω1. Now, from ( 67 ) and since r^nYnmθ,ϕ and r^-n+1Ynmθ,ϕ are orthogonal in L2Ω11, it follows

On the other hand, we have

So

Thus, from ( 98 ), ( 99 ) and ( 101 ) we get

also, from ( 94 ) we get

Making the computations, from ( 102 ) and ( 103 ) we get that anm and bnm are equivalent to n32R2R1nVnm and nR0R1nR2R1nVnm, respectively, and then anm2+bnm2 is equivalent to n3R2R12n1+1n2R0R12nVnm2. Consequently, ( 97 ) is necessary and sufficient for V to be reproducible by a source in H1Ω11.

However, these classes H1Ω11 are not unicity classes. In fact, if f ₁ and f ₂ are sources in this class generating the same measurement V in MH1Ω11, thus f1-f2∈KerA, so

Then, we get anm1=anm2, but nothing can be said about bnm1 and bnm2. From this it can be deduced that an unicity subclass contained in H1Ω11 is obtained considering harmonic functions given by series with respect only the first terms; i.e., those of the following form:

Starting from this development, and computing their coefficients via the equation Af = g as before, we get

From now, this subclass of harmonic sources will be denoted by

Note that sources class (107) precisely matches the set of sources that equal harmonic sources in Ω1 when restricted at Ω11, and vanish at Ω10.

The following result summarizes these results.

Theorem 5. If 0<R0<R1, then the admissible data sets MH1Ω11 and MH01Ω11 coincide with the admissible data set MH1Ω1. These sets are the potential distributions measurements V on the scalp whose Fourier coefficients satisfy

In addition, the class H01Ω11 is an unicity class, while H1Ω11 is not.

This last statement could result a little bit confusing; let us clarify this point. By making use of the unicity of extension of harmonic functions, one can identify a harmonic source in H01Ω11 with its harmonic extension to Ω1. In this sense, the classes H1Ω1 and H01Ω11 coincide as sets, while that H1Ω11 strictly contains H1Ω1.