Received 29 October 2014; revised 20 November 2014; accepted 8 December 2014

ABSTRACT

The inverse problem of magnetoencephalography (MEG) seeks the neuronal current within the conductive brain that generates a measured magnetic flux in the exterior of the brain-head system. This problem does not have a unique solution, and in particular, it is not even possible to identify the support of the current if it extends over a three-dimensional set. However, a localized current supported on a zero-, one- or two-dimensional set can in principle be identified. In the present work, we demonstrate an analytic algorithm that is able to recover a one-dimensional distribution of current from the knowledge of the exterior magnetic flux field. In particular, we consider a neuronal current that is supported on a small line segment of arbitrary location and orientation in space, and we reduce the identification of its characteristics to a nonlinear algebraic system. A series of numerical tests show that this system has a unique real solution. A special case is easily solved via the use of trivial algebraic operations.

Keywords:

Magnetoencephalography, Current Identification

1. Introduction

The brain is a conducting material and therefore, every generated neuronal current is accompanied by an induction current. Consequently, when we measure the magnetic flux density outside the head we actually measure the effects of both the neuronal as well as the induction current. This is the main problem with the inverse problem of magnetoencephalography, the fact that the induction current “hides” somehow the primary neuronal excitation. An excellent review of the electromagnetic activity of the human brain can be found in [1] , as well as in the book by Malmivuo and Plonsey [2] .

Exactly a hundred and sixty years ago Helmholtz [3] showed that it is not possible to recover an electric current within a conductor from knowledge of the magnetic flux generated outside the conductor. However, a complete quantitative characterization of what part of the current is possible to be identified was a topic of intense investigation during the last two decades and the main results can be found in [4] . Fokas proved that, independently of the geometry of the conductor, we cannot recover more than one out of the three functions that define the current, in the case of electroencephalography, and no more than two such functions in the case of magnetoencephalography. Even in the case that we have complete data from both modalities, still one out of the three functions is not recoverable. Another related question concerns localized neuronal currents. If the current is restricted to a small subset of the conducting brain tissue, is it possible to identify the characteristics of this current and especially its extent and its location? Albanese and Monk [5] proved that such localization is not possible. More precisely they showed that it is impossible to find the support of the current if the current occupies a three- dimensional subset of the brain. However, if the current is distributed over a surface, which is a two-dimensional subset, a curve, which is a one-dimension subset, or on isolated points, which form zero-dimensional subsets, then it is possible to identify it. It is the purpose of the present work to demonstrate that this is true for a one- dimensional current distribution. In particular, we consider a dipolar current distribution over a small line segment, and we develop an algorithm that reduces the identification of the position, the length and the orientation of the line segment, as well as the average dipolar moment of the current, to the solution of a nonlinear algebraic system. The solution of this system can be handled numerically.

2. The MEG Problem for a Single Dipole

Within the Quasi-Static Theory of Electromagnetism Magnetoencephalography [6] -[8] the magnetic field, generated by a dipolar current at the point having the moment, is given by the Geselowitz formula [9]

(1)

where u is the electric potential on the boundary S of the conducting medium representing the brain-head system. In Formula (1), denotes the exterior domain, is the constant conductivity of the brain tissue, is the magnetic permeability both inside and outside and n stands for the outward unit normal on the boundary S.

When is a sphere of radius a we know from the solution of the corresponding electroencephalography problem that the electric potential on the boundary of the sphere is given by [10] [11]

(2)

where stands for the normalized complex spherical harmonics

(3)

and denotes the Legendre functions of the first kind.

Inserting expression (2) in the Formula (1) and performing the indicated integration we can obtain the magnetic field outside the sphere. However, since the magnetic field in the exterior to the sphere is both solenoidal and irrotational it follows that there exists a scalar magnetic potential, which is also harmonic, such that [8]

(4)

Then, a series of calculations lead to the following expression for the magnetic potential [10] [11] ,

(5)

The above expression provides the magnetic potential in the exterior of the sphere due to a single current dipole. Therefore, it can be considered as the fundamental solution of the MEG problem for the spherical geometry [12] . Consequently, any discrete, or continuous, current distribution can be obtained through summation, or integration, respectively, of the above fundamental solution [13] .

3. The Field of a Linearly Distributed Current

We consider here the special case where the neuronal current is supported on a small segment of a smooth curve which is parametrically centered at the point. Let this curve be represented by the equation

(6)

The neuronal current is then described by the function,. Since the support curve is taken to be small we can approximate the current by the linear part of its Taylor expansion, that is

(7)

where the symbol denotes tensor product.

In particular, if the curve is a small line segment of length, centered at and oriented along the direction, that is

(8)

then representation (7) is written as

(9)

where provides an average moment, and provides an average directional derivative of the current along the direction.

Next we calculate the total potential which is generated by the approximate current (9). We recall that our ultimate goal is to invert the MEG data in order to identify the quantities, , and L, which are nine particular numbers, considering that the direction has two independent components. Therefore, we should be able to obtain these nine numbers from a few initial terms of the expansion (5).

Formula (5), for the excitation dipole, is written as

(10)

Using the standard expressions of the Legendre polynomials [14] and performing the indicated calculation we obtain the following relations, which are written in dyadic form [15] in order to isolate the factors that are going to be integrated

(11)

(12)

(13)

The symbol denotes the identity dyadic, “:” defines the double contraction [15]

(14)

and similarly the triple contraction is defined as

(15)

The exterior potential, given in (10), can be written in its Cartesian form [11] [13] as follows

(16)

where the coefficients

(17)

(18)

(19)

are homogeneous harmonic functions [13] .

In what follows we insert the expressions (8) and (9) in (17), (18) and (19) and integrate the resulting equations with respect to from to. Performing these calculations we arrive at the expressions

(20)

(21)

(22)

Finally, we replace the above expressions of the harmonic functions in the expansion (16) and obtain the Cartesian representation of the exterior potential up to the terms of order. That solves the relative forward MEG problem for a neuronal excitation that is supported on a small line segment.

4. Determination of the Current

The harmonic functions H₁, H₂ and H₃ are homogeneous polynomials of degrees 1, 2 and 3, respectively, that is

(23)

(24)

where, because of harmonicity, we should have the constrain

(25)

and

(26)

together with the constrains

(27)

(28)

(29)

In the idealized case where the exterior magnetic potential is known, the expansion (16) is known and therefore the coefficients, and are also known. Hence, if we rewrite the polynomials, and in terms of the Cartesian monomials that appear in (23), (24) and (26), then we can utilize their linear independence to equate each monomial with the corresponding known coefficients, or.

Equations (20) and (23) imply immediately that

(30)

Then, from Equations (30) and (33) we obtain the six relations

(31)

(32)

(33)

(34)

(35)

(36)

where it is easily shown that condition (25) holds.

Similarly, from Equations (22) and (26) we obtain

(37)

(38)

(39)

for the cubic terms, and, respectively. For the cross-terms and we obtain the expressions

(40)

and

(41)

Similarly, for the cross-terms and we obtain

(42)

and

(43)

while, for the cross-terms and we obtain

(44)

and

(45)

Finally for the product term we obtain

(46)

It is straightforward to verify that the three constrains (27)-(29) are satisfied.

The set of the 16 equations, which are the 20 scalar equations appearing in (30)-(46) minus the four constrains (25) and (27)-(29), defines a nonlinear system for the determination of the 12 independent variables, , , and, three components for each one of the vectors, , , two components for the direction vector and one for the length. In fact, we can simplify this system as follows. In view of Equation (30) the three components of the exterior product provide the relations

(47)

(48)

(49)

and these relations reduce the Equations (31)-(36) to

(50)

(51)

(52)

(53)

(54)

(55)

Furthermore, utilizing the Equations (50)-(52) we arrive at the relations

(56)

(57)

(58)

which allow rewriting Equations (37)-(46) as follows

(59)

(60)

(61)

(62)

(63)

(64)

(65)

(66)

(67)

and

(68)

Because of the constrains (27)-(29), only 7 out of the 10 equations (59)-(68) are independent. Then, the reduced set of these 7 independent equations, plus the 6 equations (47)-(49) and (53)-(55) provides a nonlinear system for the determination of the unknown quantities, , , and. However, since some of these quantities enter the system through the components of exterior products, it follows that the above quantities cannot be completely specified. For example, from Equations (47)-(49) it follows that it is not possible to identify the three components of from the exterior product of with the position vector, since the component of that is parallel to gives a vanishing term. Hence, this component of forms a “silent” source [8] . The solution of this system can easily be obtained with the use of classical computational methods.

To illustrate the inversion algorithm we consider the following special case.

Special Case. Let us assume that we have the a-priori information that the line segment is oriented along the -axis and that its middle point is. This choice leads to

(69)

(70)

(71)

(72)

(73)

(74)

Inserting the expression of in the equations for and we obtain the following system for the determination of and

(75)

(76)

which immediately gives the values of and

(77)

(78)

Then from (69) we obtain

(79)

and from (70) we obtain

(80)

Finally, from (72) and (73) we obtain a system for the unknowns and, from which we obtain

(81)

(82)

with and given by (77) and (78), respectively. Note that the only constants that remain unspecified are and, but the constant is not needed, and the constant can not be determined since the component is parallel to the position vector, and therefore their exterior product vanishes.

Acknowledgements

The present work is part of the project “Functional Brain”, which is implemented within the “ARISTEIA” Action of the “OPERATIONAL PROGRAMME EDUCATION AND LIFELONG LEARNING” and is co-funded by the European Social Fund (ESF) and National Resources.

References