Algebraic Reconstruction of Current Dipoles and Quadrupoles in Three-Dimensional Space

This paper presents an algebraicmethod for an inverse source problem for the Poisson equation where the source consists of dipoles and quadrupoles.This sourcemodel is significant in themagnetoencephalography inverse problem.Theproposedmethod identifies the source parameters directly and algebraically using data without requiring an initial parameter estimate or iterative computation of the forward solution. The obtained parameters could be used for the initial solution in an optimization-based algorithm for further refinement.


Introduction
Inverse source problems for the Poisson equation have a wide variety of applications such as bioelectromagnetic inverse problems.Magnetoencephalography (MEG) and electroencephalography (EEG) are typical examples in which the current source inside the human brain is inversely reconstructed from measurements of the magnetic field outside the head and the electric potential on the head surface [1][2][3].Due to the quasistatic property of the magnetic and electric fields, MEG and EEG exhibit high temporal resolution compared to other tools for the noninvasive inspection of brains.Hence, the development of a stable inverse algorithm which reconstructs the source with high spatial resolution is of importance for MEG and EEG.
Conventional algorithms for the MEG inverse problem, whose solution is generally nonunique [4], can be divided into two categories [2]: the parametric and imaging approaches.The parametric approach assumes that the current source inside the brain is expressed by a relatively small number of equivalent current dipoles (ECDs), which guarantees the uniqueness of the solution [4] and finds their number, positions, and dipole moments.A typical method uses optimization-based algorithms to minimize the squared sum of the difference between the forward solution and data measured at a finite number of sensor positions (see, e.g., [5]).The forward solution is given by either the Geselowitz formula, which assumes that the head consists of arbitrarliy-shaped layered domains, or the Sarvas formula, which assumes that the head consists of concentric spheres [1].A problem with optimization-based algorithms is that they require an initial solution close to the true solution, or the algorithm will often converge to a local minimum.To resolve this problem, algebraic methods have been proposed and are gaining increased importance [6][7][8][9][10].In these methods, the ECD parameters are reconstructed directly and algebraically from the boundary measurements of the magnetic field.
In contrast, the second category of inverse algorithms, the imaging approach, assumes that the many dipoles are set at the vertices of a mesh overlaying the cortical surface and solves for their dipole moments.An advantage of this approach is that the spatially distributed source can be reconstructed by solving linear equations for the unknown dipole moments.However, it has the disadvantages that the solution is not unique and an oversmoothed solution is often obtained by regularization, such as the minimum  2 norm solution [11].A method using the minimum  1 norm solution has also been proposed to reconstruct a sparse source [12].
Recently, a new method, the multipolar representation of the source, has been developed that incorporates some of the advantages of the above two methods, and has attracted considerable attention [13][14][15][16][17][18].In this model, instead of expressing the current source by an equivalent current dipole, an equivalent dipole and quadrupole [15,16] or an equivalent dipole and octopole [13,14] are used where the quadrupole or octopole is determined depending on the spatial extent of the support of the current source.Jerbi et al. showed that the centroid of the spatially distributed source, which they called a patch source, can be estimated more accurately using the dipole and quadrupole model than using the dipole model by means of the nonlinear least squares method [16].We considered a two-dimensional (2D) problem using a complex analysis framework and proposed an algebraic method to reconstruct the dipole and quadrupole parameters directly from the boundary measurements of the meromorphic function [17].The aim of the current paper is to extend our algebraic method from the 2D case to the 3D case so that the dipole and quadrupoles, which equivalently represent the neural current, can be reconstructed from the magnetic field data without an initial parameter estimate or iterative computing of the forward solution.In [19,20], we proposed an algebraic method when the dipoles were distributed in a plane parallel to the -plane, which is a very special case and is severely restricted in its practical usage.This paper derives a method for the general case.From a practical viewpoint, our method can provide an estimate of the number of patches as well as a good initial solution close to the true solution for optimization-based algorithms.
This paper is organized as follows.In Section 2, the forward problem with the dipole-quadrupole source model is summarized and the inverse problem is formulated.In Section 3, we propose a method to reconstruct the coordinates of the dipole-quadrupole source by solving simultaneous algebraic equations of second degree.A method to estimate the -coordinates is also proposed.Section 4 is devoted to numerical simulations for estimating spatially distributed dipoles.

Forward and Inverse Problems with the Dipole and Quadrupole Source Model
Let Ω 1 , Ω 2 , Ω 3 , and Ω 4 be concentric balls centered at the origin in 3D space, where Ω  ⊂ Ω +1 for  = 1, 2, 3, as shown in Figure 1.Here, Ω 1 , Ω 2 /Ω 1 , and Ω 3 /Ω 2 represent the brain, skull, and scalp, respectively.Ω 3 represents the head.We assume that the radial component of the magnetic field is measured on the sphere Γ = Ω 4 with the radius of .
Although we use this simple head model as well as the spherical sensor surface in this paper, our method can be extended to a more realistic case when the head is modeled by a piecewise homogeneous layered domain and the sensors are set on an arbitrarily shaped surface as assumed in the algebraic method for the dipole source model [10].Let J  (r) represent the neural currents in the innermost domain Ω 1 .The magnetic permeability is assumed to be constant  0 in the whole space.First, we derive the representation of the radial magnetic field in terms of the equivalent dipoles and quadrupoles, which has already been given in [15], to determine how the quadrupole terms equivalently express the spatial extent of the source.
It is known that with this head model, the radial component of the magnetic field at r is equal to that when J  (r) exists in infinite homogeneous space and is given by the Biot-Savart law: Let us assume that the neural activity is localized at several domains   (⊂ Ω 1 ) for  = 1, 2, . . .,  so that J  (r) is given by where    is the characteristic function.Then, the radial component of the magnetic field is given by Assume that   is contained in a ball   ⊂ Ω 1 , and there exists a "centroid" r  of   such that From the Taylor series expansion, we have ( Substituting this expansion into the right-hand side of (3) gives where   is the cross-product tensor defined by and hence is written as In (6), ":" represents the tensor contraction defined by  :  = ∑ ) , then we have Note here that we define   to be symmetric,  , =  , , since in (6) only the symmetric part of ∫   j  (r  )(r  − r  )  contributes to r ⋅ B as shown in [15].Equation (10) was given in (46) in [15] (when  = 1).The first term in (10) is the magnetic field created by an "equivalent current dipole" p  , and the second term is the magnetic field created by an "equivalent current quadrupole"   .Note here that p  does not depend on the size of   , while   depends on the spatial extent of j  in   .  is a 3 × 3 tensor of order 2 and is called the quadrupole moment tensor.
Here, we can show that the truncation of (10) up to the second term, is the forward solution for the source model given by where  is the Dirac delta function (see Appendix B).We call (12) the "equivalent dipole and quadrupole source model" or simply the "dipole-quadrupole source model." We are now ready to state our inverse problem.
Problem 1. Assume that a head Ω 3 is modeled by a spherically layered domain.Assume that the neural current is expressed by (12) supported in the inner-most ball Ω 1 representing the brain.Then, reconstruct the number  and positions r  of the dipole-quadrupole source in (12) from the radial component of the magnetic field measured on the spherical boundary Γ which encloses Ω 3 .
Remark 1. Once the number and positions of the dipolequadrupole source are determined, the dipole moments p  and quadrupole moment tensors   can be obtained by solving the relevant linear equation (11).In this paper, we concentrate on reconstructing  and r  .

Algebraic Reconstruction Method
According to (88) in [15], we can express (1) by the multipole expansion where and    (cos ) are the associated Legendre polynomials.As shown in (84) in [15], the multipole moment has the following relationship with J  : Substituting ( 12) into (15) and using (B.1) gives On the other hand, using the orthogonality of the spherical harmonics, the multipole moment is shown to have another relationship with the radial magnetic field on a spherical surface Γ [21]: where n  is the outward unit normal vector at r  on Γ. Equating ( 16) and ( 17) for  ≥  ≥ 0 gives algebraic equations relating the unknown parameters of the sources to the radial magnetic field on Γ.Although (17) requires the radial components of the magnetic field on the whole boundary Γ, Taulu et al. proposed a new approach for obtaining   from the data on the part of Γ by using the truncated spherical harmonic expansion [22].Using the method, we have algebraic equations about the unknown source parameters with the practically measurable MEG data.

Algebraic Reconstruction of the 𝑥and 𝑦-Coordinates.
As we derived the algebraic method for the dipole source model in [10], we now use the multipole moments of the sectorial harmonics represented by  =  = , where  ∈ N. From the relationship     (, ) = (2 − 1)!!( + )  [23], we have The primes in both sides are removed for simplicity.Since it holds that we obtain where See Appendix C for the derivation of (21).Substituting ( 20) and ( 21) into (18) gives We now put and rewrite (22) as Comparing with the special case where the dipoles were distributed in a plane parallel to the -plane in [20], one finds that   has an extra term: ( , −  , )(  +   ).Equation ( 24) has the same form as that of (6) in [17], and consequently the dipole-quadrupole position projected on the -plane,   ( = 1, 2, . . ., ), can be algebraically reconstructed from   ( = 0, 1, . . ., 2 − 1), which is shown as follows.In Section 3.1, we assume that  is known a priori and construct an algorithm to estimate   in (24).A method to estimate  as well as   and   is then described in Section 3.2.
By definition,  (2)  , are the Hankel matrices expressed by Note that the components along the first and second lines in the antidiagonal direction become zero.Repeating this procedure to eliminate the components on the th line in the anti-diagonal direction for  = 3, . . ., 2 − 1, we obtain where  () , are the Hankel matrices whose components along the first through th lines in the anti-diagonal direction are zero.Here, (34) for  = 2 − 1 is a linear equation for   .Solving it and substituting the solution into (34) for  = 2−2 gives a linear equation for  −1 .Repeating the backward substitution gives  −2 , . . .,  1 , successively.The conditions for the unique solvability of these linear equations are   ̸ =   for  ̸ =  and   ̸ = 0. See [17] for more details.Once  1 , . . .,   are obtained, we have   as the roots of   + 1  −1 +⋅ ⋅ ⋅+  = 0.
In a practical situation when data includes noise, the above procedure proposed in [17] to transform the seconddegree equations (26) into linear equations (34) to obtain   , . . .,  1 may be sensitive to the noise due to the cancellation.To reduce the cancellation error, we propose the following algorithm.First, we solve the second-degree equations (26) for  = 0, . . .,  − 1 using a well-established method for solving simultaneous algebraic equations, for example, by means of the Gr ö bner bases.Then, out of the obtained 2  solutions, ( 1, , . . .,  , ) for  = 1, . . ., 2  , we choose one which minimizes the sum of the absolute values of the lefthand sides of ( 26) for  = , . . ., 2 − 1 given by In other words, to obtain the theoretically unique solution of ( 26) for  = 0, . . ., 2 − 1, first we find 2  candidates from (26) for  = 0, . . .,  − 1 and then choose the one which best explains the remaining equations ( 26) for  = , . . ., 2 − 1.

3.2.
Reconstruction of   ,   , and .To estimate , following the method for the dipole source model [10], we assume that there are   (>N) dipole-quadrupole sources and then estimate   for  = 1, 2, . . .,   using the algorithm in Section 3.1.
Once they are obtained,   and   for  = 1, 2, . . .,   are linearly solved using (24) for  = 0, 1, . . ., 2  − 1.Then, we compute | +1 /  | and | +1 /  | for  = 1, 2, . . .,   −1, which are expected to be sufficiently small when  = .Practically, we estimate  such that these ratios become smaller than some threshold set a priori.The thresholds should be determined by the ratios of the noise level contained in the data to the dipole and quadrupole strength which can be regarded as a physiologically meaningful source.As for the dipole source model in the 2D problem, the threshold is theoretically evaluated in the context of the Padé approximation [24].
A similar theory for the dipole-quadrupole source model, although greatly required, is left for further research; in this paper we show only numerical simulations in Section 4.

Numerical Simulations
In this section, our algorithm is verified numerically.Γ was assumed to be a sphere with  = 0.12 m on which 361 measurement points were distributed uniformly using the spherical t-design [25].To validate our algorithm for the dipolequadrupole model, we assumed that the data was available on the whole sphere which enclosed the source.Identification using data on a part of Γ is an important aspect of future studies, for which the method proposed by Taulu et al. [22] would be useful.

Demonstration of Our Method.
In this subsection, We examine the case where there are  = 2 dipole-quadrupole sources whose parameters are given in Table 1.The theoretical data was computed using (11) to which Gaussian noise with  = 10 −2 was added where the noise level  is defined by the ratio of the standard deviation of the noise to the root mean squares of the data.100 data sets with the different noise added were used for reconstruction.
First, to determine the number of the sources , we assume that there are   = 3 dipole-quadrupole sources.
Mathematical Problems in Engineering  Figure 3 shows the localization result projected on the and -planes.We find that, under this noise level, two dipole-quadrupole sources were stably reconstructed.The average and standard deviations of the computational time to identify the parameters of the two dipole-quadrupole sources were 0.69 sec and 0.01 sec, respectively, using an Intel Core i7 CPU 870 (2.93 GHz) with 8 GB RAM, showing that our algebraic algorithm can determine the source parameters with low computational cost.

Effect of Noise on Localization Accuracy.
To examine the stability of our algorithm, the noise level  was varied in the wide range 10 −8 ≤  ≤ 10 −1 .The source and sensor configuration were the same as in Section 4.1.Figure 4 shows    the standard deviation (s.d.) of the localization errors with respect to the noise level  for the two poles.It is observed that s.d.increases as  increases.For this source in the sphere with the radius 0.12 m, to realize s.d.less than 0.01 m, one finds that the noise level should be less than  = 10 −2 .

Effect of the Distance between Two Sources on Localization
Accuracy.We varied the distance between the two dipolequadrupole sources: the two dipole-quadrupole source were set at (0.02, ±, 0.07) where 2 was changed from 0.01 through 0.09. Figure 5 shows the relative localization error (error divided by the radius of Γ, 0.12).It is observed that the distance between the two dipole-quadrupole sources should be larger than about 0.06, that is the half of the radius of Γ, in order to obtain the relative error less than 10 −2 .

Example
When  = 3.We examine the case when there are  = 3 dipole-quadrupole sources whose parameters are given in Table 2.The noise level was set at  = 10 −2 .In this subsection, we assumed that  is known a priori.Figure 6 shows the localization result.Comparing with the result when  = 2 in Figure 3, the biases as well as the standard deviations of the localization errors become  larger when  = 3.It is numerically shown in [26] that, when using dipole-only source model, stability gets worse when  increases and it is difficult to estimate the dipoles greater than about  = 5.The example in this subsection shows that it is more difficult to estimate the large number ( ≥ 3) of dipole-quadrupole sources.This is considered due to fast spatial decay of the quadrupolar field.However, the dipole-quadrupole source even with  = 1 or 2 has an advantage that it can well estimate a spatially distributed sources, especially when they are modeled with two, close, oppositely directed dipoles, which is a weak point of the conventional dipole model.This is shown in the next subsection.

Identification of Dipoles Distributed on Curved Surfaces.
In this subsection, we compare our algebraic method assuming the dipole-quadrupole model (DQM) with a conventional algebraic method assuming the dipole model (DM) for estimating the spatially distributed dipole sources.To model dipoles on cerebral convolutions, we assume that dipoles are placed on a mesh on a half-cylinder with a radius of  = 5 mm and a height of ℎ = 5mm, as shown in Figure 7.There are six dipoles in the circumferential direction by five in   the longitudinal direction, and hence a total of 30 dipoles on the half-cylinder.All the dipoles are aligned perpendicular to the surface of the cylinder to model the fact that the dipoles are perpendicular to the cerebral surface.We examined the following two cases.
Case 1.A single half-cylinder source at r 1 = (20, 50, 20) mm.The vectors which determine the posture of the cylinder, e 12 and e 11 , are set to be e  and e  , respectively, where e  and e  are the unit vectors in the  and  directions at r 1 .(see Figure 8(a)).In this case, the total dipole moment p 1 is almost parallel to r 1 ; the angle between them is 2.4 degrees.Since a radial dipole is a silent source for the radial component of the magnetic field [1], this cylindrical source is regarded as being almost quadrupolar.
Case 2. The same half-cylinder at r 1 as in Case 1, but with e 12 = −e  and e 11 = e  (Figure 8(b)).In this case, the angle between r 1 and p 1 is 87 degrees, so that the source has a detectable equivalent dipole moment as well as the equivalent quadrupole moment.
We computed the forward solution generated by the 30 elemental dipoles using (1).Note that (11) was not used to compute the theoretical data.  = 3 was assumed to estimate .

Case 1.
Figure 9 shows the estimated |  | normalized by max{  } using DM.From the fact that | 3 / 2 | = 0.005,  is estimated to be two using DM.The reason why not a single equivalent dipole but two dipoles are estimated for a single cylindrical source is that the sum of the dipoles on the half-cylindrical surface is directed to the radial direction which is "silent" to the radial magnetic field outside the head.In contrast, Figures 10(a Figure 11 shows the localization result using DM and DQM.It is observed that the two positions estimated with DM were far apart from the cylindrical surface, whereas DQM well estimated the center of the cylinder.This is a great   advantage of DQM; although it is difficult for DM to estimate stably two, close, oppositely directed dipoles which are tangential to the spherical surface, it is easy for DQM.  13 using DQM, both model estimate  = 1.This is because the total of the elemental dipoles on the cylinder is almost perpendicular to the radial direction in Case 2 so that this cylindrical source can be regarded as a single equivalent current dipole.
Figure 14 shows the localization result.Both DM and DQM can well estimate the center position of the cylindrical source in this case.

Conclusion
We considered an inverse source problem of the Poisson equation for the radial component of the magnetic flux density when the current source is equivalently represented by multiple dipoles and quadrupoles in layered, spherical domains.By expressing the magnetic field in terms of the spherical harmonic expansion, we showed that the sectorial harmonics gave the simultaneous algebraic equations relating the dipole and quadrupole parameters to the boundary data.The equations had the same form as derived in the simple-and double-pole reconstruction problem in the 2D case [17], so that the source parameters can be algebraically reconstructed.We proposed a method to obtain the and -coordinates of the dipole-quadrupole source by solving the simultaneous second-degree equations, with the result of which the -coordinates are also determined.Numerical simulations show that the proposed algorithm estimates the number and centroid positions of the spatially extended dipoles well, especially when they include elemental dipoles oriented in opposite directions and their equivalent dipole moment is almost parallel to the radial direction of the spherical domain.A theoretical analysis to determine the thresholds depending on the noise in the algorithm for the estimation of the number of the dipole-quadrupole sources as well as an extension of the method to the case when the data is available only for part of the spherical boundary is left for further research.
C. Derivation of (21) It holds that

Figure 2
Figure 2 shows the obtained |  | and |  | for  = 1, 2, 3.One observes that | 3 | and | 3 | are much smaller than |  | and |  | for  = 1, 2, respectively.In fact, the geometric means of | 3 / 2 | and | 3 / 2 | for 100 data sets were 3.8 × 10 −3 and 7.2 × 10 −3 , respectively.From this, we can judge that there are substantially  = 2 poles.Figure3shows the localization result projected on the and -planes.We find that, under this noise level, two dipole-quadrupole sources were stably reconstructed.The average and standard deviations of the computational time to

Figure 4 :
Figure 4: Standard deviation of the localization errors with respect to the noise level .

Figure 5 :
Figure 5: Relative localization error with respect to the distance between the two dipole-quadrupole sources.

Figure 11 :
Figure 11: Localization result of the cylindrically distributed dipole sources using dipole model (DM) and dipole-quadrupole model (DQM).