On the Inverse MEG Problem with a 1-D Current Distribution

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-, oneor 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.


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 com-plete 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 threedimensional 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 onedimensional 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.

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 0 r having the moment Q , is given by the Geselowitz formula [9] ( ) ( ) ( ) ( ) where u is the electric potential on the boundary S of the conducting medium Ω representing the brain-head system.In Formula (1), C Ω denotes the exterior domain, σ is the constant conductivity of the brain tissue, 0 µ 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] ( ) where m n Y stands for the normalized complex spherical harmonics and m n P 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 B in the exterior to the sphere is both sole- noidal and irrotational it follows that there exists a scalar magnetic potential U , which is also harmonic, such that [8] ( ) ( ) Then, a series of calculations lead to the following expression for the magnetic potential [10] [11], ( ) ( ) The above expression provides the magnetic potential in the exterior of the sphere due to a single current dipole { } 0 , r Q .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].

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 0 r .Let this curve be represented by the equation 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 where the symbol ⊗ denotes tensor product.
In particular, if the curve is a small line segment of length 2L , centered at 0 r and oriented along the direction ( ) then representation ( 7) is written as , , Q J r 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 Q , 0 r , α and L, which are nine par- ticular 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 r r J r r r r r r r r (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 ( ) ( ) ) ( ) The symbol  I denotes the identity dyadic, ":" defines the double contraction [15] ( ) ( ) ( )( ) (14) and similarly the triple contraction is defined as The exterior potential, given in (10), can be written in its Cartesian form [11] [13] as follows ; 4π where the coefficients 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 t from L − to L .Performing these calculations we arrive at the expressions Finally, we replace the above expressions of the harmonic functions 1 2 3 , , H H H in the expansion (16) and obtain the Cartesian representation of the exterior potential U up to the terms of order 5  r − .That solves the relative forward MEG problem for a neuronal excitation that is supported on a small line segment.

Determination of the Current
The harmonic functions H 1 , H 2 and H 3 are homogeneous polynomials of degrees 1, 2 and 3, respectively, that is ( ) where, because of harmonicity, we should have the constrain together with the constrains In the idealized case where the exterior magnetic potential U is known, the expansion ( 16) is known and therefore the coefficients A , B and C are also known.Hence, if we rewrite the polynomials 1 H ,

2
H and 3 H in terms of the Cartesian monomials that appear in (23), ( 24) and ( 26), then we can utilize their linear in- dependence to equate each monomial with the corresponding known coefficients A , B or C .Equations ( 20) and ( 23) imply immediately that ( ) ( ) ( ) Then, from Equations ( 30) and ( 33) we obtain the six relations where it is easily shown that condition (25) holds.Similarly, from Equations ( 22) and (26) we obtain for the cubic terms 3 1 x , 3   2   x and 3   3   x , respectively.For the cross-terms 2 1 2 x x and 2 1 2 x x we obtain the expressions Similarly, for the cross-terms 2 2 3 x x and 2 2 3 x x we obtain while, for the cross-terms 2  ) ( Finally for the product term 1 2 3 x x x we obtain ( )( ) 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 0 r , â , l , Q and L , three components for each one of the vectors 0 r , l , Q , two components for the direction vector â and one for the length L .In fact, we can simplify this system as follows.In view of Equation (30) the three components of the exterior product 0 × Q r provide the relations ( ) Furthermore, utilizing the Equations ( 50)-( 52) we arrive at the relations ( ) ( ) ( ) which allow rewriting Equations ( 37)-(46) as follows 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 0 r , â , Q , l and L .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 Q from the exterior product of Q with the position vector 0 r , since the component of Q that is parallel to 0 r gives a vanishing term.Hence, this component of Q 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 1 x -axis and that its middle point is ( ) . This choice leads to Finally, from (72) and ( 73) we obtain a 2 2 × system for the unknowns 1 Q and 3 l , from which we obtain