#### Abstract

This paper presents an algebraic method for an inverse source problem for the Poisson equation where the source consists of dipoles and quadrupoles. This source model is significant in the magnetoencephalography inverse problem. The proposed method 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.

#### 1. 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–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–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 norm solution [11]. A method using the minimum 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–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.

#### 2. Forward and Inverse Problems with the Dipole and Quadrupole Source Model

Let , , , and be concentric balls centered at the origin in 3D space, where for , as shown in Figure 1. Here, , , and represent the brain, skull, and scalp, respectively. represents the head. We assume that the radial component of the magnetic field is measured on the sphere 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 represent the neural currents in the innermost domain . The magnetic permeability is assumed to be constant 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 is equal to that when 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 for so that 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 , and there exists a “centroid” 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 for second-order tensors and whose components are given by and , respectively. See Appendix A for the derivation of (6).

Now, let then we have Note here that we define to be symmetric, , since in (6) only the symmetric part of contributes to as shown in [15]. Equation (10) was given in (46) in [15] (when ). The first term in (10) is the magnetic field created by an “equivalent current dipole” , and the second term is the magnetic field created by an “equivalent current quadrupole” . Note here that does not depend on the size of , while depends on the spatial extent of in . is a 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 is modeled by a spherically layered domain. Assume that the neural current is expressed by (12) supported in the inner-most ball representing the brain. Then, reconstruct the number and positions of the dipole-quadrupole source in (12) from the radial component of the magnetic field measured on the spherical boundary which encloses .

*Remark 1. *Once the number and positions of the dipole-quadrupole source are determined, the dipole moments and quadrupole moment tensors can be obtained by solving the relevant linear equation (11). In this paper, we concentrate on reconstructing and .

#### 3. Algebraic Reconstruction Method

According to (88) in [15], we can express (1) by the multipole expansion where and are the associated Legendre polynomials. As shown in (84) in [15], the multipole moment has the following relationship with : Substituting (12) into (15) and using 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 is the outward unit normal vector at on . Equating (16) and (17) for 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.

##### 3.1. 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 . From the relationship [23], we have The primes in both sides are removed for simplicity. Since it holds thatwe obtain where represents -component of + -component of . 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 in [17], and consequently the dipole-quadrupole position projected on the -plane, , can be algebraically reconstructed from , 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.

Let us define by the coefficients of It was proved in [17] that (24) can be transformed into the second-degree equations for : where and is the Hankel matrix given by Moreover, it has been shown that the equations (26) for can be turned into linear equations as follows. First, from those equations, we make equations where By this definition, are the Hankel matrices whose -components are zero. We denote their components by Next, from those equations, we make equations: where By definition, 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 , 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 is a linear equation for . Solving it and substituting the solution into (34) for gives a linear equation for . Repeating the backward substitution gives , successively. The conditions for the unique solvability of these linear equations are for and . See [17] for more details. Once are obtained, we have as the roots of .

In a practical situation when data includes noise, the above procedure proposed in [17] to transform the second-degree equations (26) into linear equations (34) to obtain 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 using a well-established method for solving simultaneous algebraic equations, for example, by means of the Gröbner bases. Then, out of the obtained solutions, for , we choose one which minimizes the sum of the absolute values of the left-hand sides of (26) for given by In other words, to obtain the theoretically unique solution of (26) for , first we find candidates from (26) for and then choose the one which best explains the remaining equations (26) for .

##### 3.2. Reconstruction of , , and

To estimate , following the method for the dipole source model [10], we assume that there are dipole-quadrupole sources and then estimate for using the algorithm in Section 3.1. Once they are obtained, and for are linearly solved using (24) for . Then, we compute and for , 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.

##### 3.3. Algebraic Reconstruction of the -Coordinates

To determine the -coordinates of the dipole-quadrupole positions, after we determine the -coordinates of the source, we use the method proposed in [7]. We put and in (16) and (17), where . From the identity [23], we have Now, let Then, from × (18) and × (36), we obtain respectively. In the left-hand sides of (38) and (39), since from (37) it holds that for all irrespective of whether or not is equal to . Also, we have where is the Kronecker delta, and See Appendix D for the derivation of (42) and (43).

Substituting (40) through (43) into (38) and (39) gives from which we obtain This is the reconstruction formula of the -coordinates: once for are obtained, then are determined using (46) for each with defined by (37).

##### 3.4. Algorithm

Our algorithm is summarized as follows.

*Step 1 (estimate the number of the dipole-quadrupole sources). *(1) Assume that there are dipole-quadrupole sources. Using the algorithm in Step 2 (1)–(5) below, where is replaced by , obtain and for .

(2) Compute and for . When and for the thresholds and , which are appropriately set *a priori*, we estimate to be .

*Step 2 (estimate the position parameters of dipole-quadrupole sources). *(1) Given on , compute in (23) for and construct the Hankel matrices (27) for .

(2) Solve the simultaneous second-degree equation (26) for , using the Gröbner bases to get solutions.

(3) Out of the solutions, choose one solution that minimizes (35).

(4) Solve to obtain .

(5) Solve (24) for and , where .

(6) Using the obtained , compute in (37) for and determine using (46).

#### 4. Numerical Simulations

In this section, our algorithm is verified numerically. was assumed to be a sphere with m on which 361 measurement points were distributed uniformly using the spherical t-design [25]. To validate our algorithm for the dipole-quadrupole 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.

##### 4.1. Demonstration of Our Method

In this subsection, We examine the case where there are dipole-quadrupole sources whose parameters are given in Table 1. The theoretical data was computed using (11) to which Gaussian noise with 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 dipole-quadrupole sources. Figure 2 shows the obtained and for . One observes that and are much smaller than and for , respectively. In fact, the geometric means of and for 100 data sets were and , respectively. From this, we can judge that there are substantially poles.

**(a)**

**(b)**

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.

**(a)**

**(b)**

##### 4.2. Effect of Noise on Localization Accuracy

To examine the stability of our algorithm, the noise level was varied in the wide range . 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 .

##### 4.3. Effect of the Distance between Two Sources on Localization Accuracy

We varied the distance between the two dipole-quadrupole sources: the two dipole-quadrupole source were set at where was changed from through . 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 .

##### 4.4. Example When

We examine the case when there are dipole-quadrupole sources whose parameters are given in Table 2. The noise level was set at . In this subsection, we assumed that is known *a priori*. Figure 6 shows the localization result. Comparing with the result when in Figure 3, the biases as well as the standard deviations of the localization errors become larger when . 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 . The example in this subsection shows that it is more difficult to estimate the large number of dipole-quadrupole sources. This is considered due to fast spatial decay of the quadrupolar field. However, the dipole-quadrupole source even with 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.

**(a)**

**(b)**

##### 4.5. 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 mm and a height of mm, 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 mm. The vectors which determine the posture of the cylinder, and , are set to be and , respectively, where and are the unit vectors in the and directions at . (see Figure 8(a)). In this case, the total dipole moment is almost parallel to ; 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.

**(a)**

**(b)**

*Case 2. *The same half-cylinder at as in Case 1, but with and (Figure 8(b)). In this case, the angle between and 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. was assumed to estimate .

###### 4.5.1. Case 1

Figure 9 shows the estimated normalized by using DM. From the fact that , 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) and 10(b) show the estimated normalized by and normalized by , respectively, using DQM. From and , is estimated to be one.

**(a)**

**(b)**

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.

**(a)**

**(b)**

###### 4.5.2. Case 2

Figures 12, 13(a), and 13(b) show the estimated using DM and and using DQM, respectively. In this Case, from in Figure 12 using DM and and in Figure 13 using DQM, both model estimate . 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.

**(a)**

**(b)**

**(a)**

**(b)**

#### 5. 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.

#### Appendices

#### A. Derivation of (6)

Substituting the Taylor expansion of into the right-hand side of (3), we have The first term of (A.1) is rewritten as The second term of (A.1) is rewritten as Here, for a matrix whose row vectors are , , and , it holds that Using this property of the tensor contraction, (A.3) is further rewritten as Therefore, substituting (A.2) and (A.5) into (3) gives (6).

#### B.

This can be shown as follows. It is easy to obtain the dipole terms. For the quadrupole terms, we use the identity: for an arbitrary vector field , where represents the transpose. When inserting the quadrupole terms in (12) into (1), we have from (B.1) Here, it holds that and hence Thus we have (11).

#### C. Derivation of (21)

It holds thatwhich gives (21).

#### D. Derivation of (42) and (43)

It is easy to see that where , , and represent , , and , respectively. From the definition of in (37), one easily finds that for all , and hence Also, from the definition of , we have where is Kronecker’s delta, which gives (42). In the same way, which gives (43).

#### Acknowledgments

This research has been partially supported by the Moritani Scholarship Foundation and the 2012 Inamori Research Grants Program.