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

Authors: Takaaki Nara
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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 minimumL2 norm solution [11]. A method using the minimum L1 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 xy-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 Section2, 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 xy-coordinates of the dipole-quadrupole source by solving simultaneous algebraic equations of second degree. A method to estimate the z-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Ω1, Ω2, Ω3, and Ω4 be concentric balls centered at the origin in 3D space, where Ωi⊂Ωi+1 for i=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 R. 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 Jp(r) represent the neural currents in the innermost domain Ω1. The magnetic permeability is assumed to be constant μ0 in the whole space.Figure 1
Head and sensor model.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 atr is equal to that when Jp(r) exists in infinite homogeneous space and is given by the Biot-Savart law:
(1)r·B(r)=μ04π∫Ω1‍(∇′1|r-r′|×r)·Jp(r′)dv′.
Let us assume that the neural activity is localized at several domains Dk(⊂Ω1) for k=1,2,…,N so that Jp(r) is given by
(2)Jp(r)=∑k=1N‍jp(r)χDk,
where χDk is the characteristic function. Then, the radial component of the magnetic field is given by
(3)r·B(r)=μ04π∑k=1N‍∫Dk‍(∇′1|r-r′|×r)·jp(r′)dv′.
Assume that Dk is contained in a ball Bk⊂Ω1, and there exists a “centroid” rk of Dk such that
(4)|r-rk|>|r′-rk|,∀r∈Γ,r′∈Dk.
From the Taylor series expansion, we have
(5)∇′1|r-r′|=∇′1|r-r′||r′=rk+∇′∇′1|r-r′||r′=rk(r′-rk)+⋯.
Substituting this expansion into the right-hand side of (3) gives
(6)r·B(r)=μ04π∑k=1N‍r-rk|r-rk|3·(rk×∫Dk‍jp(r′)dv′)+μ04π∑k=1N‍3(r×rk)(r-rk)+|r-rk|2Xr|r-rk|5:∫Dk‍jp(r′)(r′-rk)dv′+⋯,
where Xr is the cross-product tensor defined by
(7)Xra=r×a,
and hence is written as
(8)Xr=(0-zyz0-x-yx0)forr=(x,y,z)T.
In (6), “:” represents the tensor contraction defined by A:B=∑i=13‍∑j=13‍AijBij for second-order tensors A and B whose (i,j) components are given by Aij and Bij, respectively. See Appendix A for the derivation of (6).Now, let(9)pk≡∫Dk‍jp(r′)dv′,Qk≡12(∫Dk‍jp(r′)(r′-rk)+(r′-rk)jp(r′))dv′=(qxx,kqxy,kqxz,kqyx,kqyy,kqyz,kqzx,kqzy,kqzz,k),
then we have
(10)r·B(r)=μ04π∑k=1N‍(rk×pk)·r-rk|r-rk|3+μ04π∑k=1N‍Qk:3(r×rk)(r-rk)+|r-rk|2Xr|r-rk|5+⋯.
Note here that we define Qk to be symmetric, qij,k=qji,k, since in (6) only the symmetric part of ∫Dk‍jp(r′)(r′-rk)dv′ contributes to r·B as shown in [15]. Equation (10) was given in (46) in [15] (when N=1). The first term in (10) is the magnetic field created by an “equivalent current dipole” pk, and the second term is the magnetic field created by an “equivalent current quadrupole” Qk. Note here that pk does not depend on the size of Dk, while Qk depends on the spatial extent of jp in Dk. Qk 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,
(11)r·B(r)=μ04π∑k=1N‍((rk×pk)·(r-rk)|r-rk|3r·B(r)=μ04π∑k=1N‍+Qk:3(r×rk)(r-rk)+|r-rk|2Xr|r-rk|5),
is the forward solution for the source model given by
(12)Jp(r)=∑k=1N‍pkδ(r-rk)+∑k=1N‍Qk∇δ(r-rk),
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 N and positions rk 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 dipole-quadrupole source are determined, the dipole momentspk and quadrupole moment tensors Qk can be obtained by solving the relevant linear equation (11). In this paper, we concentrate on reconstructing N and rk.


## 3. Algebraic Reconstruction Method

According to (88) in [15], we can express (1) by the multipole expansion
(13)r·B(r)=μ0∑l=0∞‍∑m=-ll‍l+12l+1MlmY^lm*(θ,ϕ)rl+1,
where
(14)Y^lm(θ,ϕ)=2l+14π(l-m)!(l+m)!×Ylm(θ,ϕ):normalizedsphericalharmonics,Ylm(θ,ϕ)=Plm(cosθ)eimϕ:sphericalharmonics,
and Plm(cosθ) are the associated Legendre polynomials. As shown in (84) in [15], the multipole moment has the following relationship with Jp:
(15)Mlm=1l+1∫Ω1‍[∇′r′lY^lm(θ′,ϕ′)]·(r′×Jp(r′))d3r′.
Substituting (12) into (15) and using (B.1) gives
(16)Mlm=1l+1∑k=1N‍((∇′(∇′r′lY^lm(θ′,ϕ′)×r′))T|r′=rk∇′r′lY^lm(θ′,ϕ′)|r′=rk+Qk:(∇′(∇′r′lY^lm(θ′,ϕ′)×r′))T|r′=rk).
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]:
(17)Mlm=2l+1μ0(l+1)∫Γ‍n′·B(r′)RlY^lm(θ′,ϕ′)dS′,
where n′ is the outward unit normal vector at r′ on Γ. Equating (16) and (17) for l≥m≥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 Mlm 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 thex- and y-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 l=m=n, where n∈ℕ. From the relationship rnYnn(θ,ϕ)=(2n-1)!!(x+iy)n [23], we have
(18)∑k=1N‍((∇(∇(x+iy)n×r))T|r=rk∇(x+iy)n|r=rk·(rk×pk)∑k=1N‍+Qk:(∇(∇(x+iy)n×r))T|r=rk)=2n+1μ0∫Γ‍n·B(r)(x+iy)ndS.
The primes in both sides are removed for simplicity. Since it holds that(19)∇(x+iy)n=n(1i0)(x+iy)n-1,(∇(∇(x+iy)n×r))T=n(i(n-1)z(x+iy)n-2-(n-1)z(x+iy)n-2i(x+iy)n-1-(n-1)z(x+iy)n-2-i(n-1)z(x+iy)n-2-(x+iy)n-1-in(x+iy)n-1n(x+iy)n-10),we obtain
(20)∇(x+iy)|r=rk·(rk×pk)=n[rk×pk]x+iy(xk+iyk)n-1,(21)Qk:(∇(∇(x+iy)n×r))T|r=rk=n(n-1)((i(qxx,k-qyy,k)-2qxy,k)zk=n(n-1)+(qyz,k-iqxz,k)(xk+iyk))(xk+iyk)n-2,
where [rk×pk]x+iy represents (x-component of rk×pk) + i×(y-component of rk×pk). See Appendix C for the derivation of (21). Substituting (20) and (21) into (18) gives
(22)∑k=1N‍(+(qyz,k-iqxz,k)(xk+iyk))(xk+iyk)n-2[rk×pk]x+iy(xk+iyk)n-1+(n-1)×((i(qxx,k-qyy,k)-2qxy,k)zk+(qyz,k-iqxz,k)(xk+iyk))(xk+iyk)n-2)=2n+1μ0n∫Γ‍n·B(x+iy)ndS.
We now put
(23)ζk=xk+iyk,μk=[rk×pk]x+iy,νk=(i(qxx,k-qyy,k)-2qxy,k)zk+(qyz,k-iqxz,k)(xk+iyk),cn=2n+3μ0(n+1)∫Γ‍n·B(x+iy)n+1dS
and rewrite (22) as
(24)∑k=1N‍μkζkn+n∑k=1N‍νkζkn-1=cn(n∈ℕ*).
Comparing with the special case where the dipoles were distributed in a plane parallel to the xy-plane in [20], one finds that νk has an extra term: (qyz,k-iqxz,k)(xk+iyk).Equation (24) has the same form as that of (6) in [17], and consequently the dipole-quadrupole position projected on the xy-plane, ζk(k=1,2,…,N), can be algebraically reconstructed from cn(n=0,1,…,2N-1), which is shown as follows. In Section 3.1, we assume that N is known a priori and construct an algorithm to estimate ζk in (24). A method to estimate N as well as μk and νk is then described in Section 3.2.Let us defineσ1,σ2,…,σN by the coefficients of
(25)∏k=1N‍(ζ-ζk)=ζN+σ1ζN-1+⋯+σN.
It was proved in [17] that (24) can be transformed into the second-degree equations for σ1,σ2,…,σN:
(26)σTHN,mσ=0(m∈ℕ*),
where σ=(σN,…,σ1,1)T and HN,m is the (N+1)×(N+1) Hankel matrix given by
(27)HN,m=(cmcm+1⋯cm+Ncm+1cm+2⋯cm+N+1⋮⋮⋮⋮cm+Ncm+N+1⋯cm+2N).
Moreover, it has been shown that the 2N equations (26) for m=0,1,…,2N-1 can be turned into linear equations as follows. First, from those 2N equations, we make 2N-1 equations
(28)σTHN,m(1)σ=0(m=0,1,…,2N-2),
where
(29)HN,m(1)≡cm+1HN,m-cmHN,m+1.
By this definition, HN,m(1) are the Hankel matrices whose (1,1)-components are zero. We denote their components by
(30)HN,m(1)=(0cm+1(1)⋯cm+N(1)cm+1(1)cm+2(1)⋯cm+N+1(1)⋮⋮⋮⋮cm+N(1)cm+N+1(1)⋯cm+2N(1)).
Next, from those 2N-1 equations, we make 2N-2 equations:
(31)σTHN,m(2)σ=0(m=0,1,…,2N-3),
where
(32)HN,m(2)=cm+2(1)HN,m(1)-cm+1(1)HN,m+1(1).
By definition, HN,m(2) are the Hankel matrices expressed by
(33)HN,m(2)=(00cm+2(2)⋯cm+N(2)0cm+2(2)⋯⋯cm+N+1(2)cm+2(2)⋯⋯⋯cm+N+1(2)⋮⋮⋮⋮⋮cm+N(2)cm+N+1(2)⋯⋯cm+2N(2)).
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 jth line in the anti-diagonal direction for j=3,…,2N-1, we obtain
(34)σTHN,m(j)σ=0(m=0,1,…,2N-1-j),
where HN,m(j) are the Hankel matrices whose components along the first through jth lines in the anti-diagonal direction are zero. Here, (34) for j=2N-1 is a linear equation for σN. Solving it and substituting the solution into (34) for j=2N-2 gives a linear equation for σN-1. Repeating the backward substitution gives σN-2,…,σ1, successively. The conditions for the unique solvability of these linear equations are ζi≠ζj for i≠j and νk≠0. See [17] for more details. Once σ1,…,σN are obtained, we have ζk as the roots of ζN+σ1ζN-1+⋯+σN=0.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 σN,…,σ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 m=0,…,N-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 2N solutions, (σ1,j,…,σN,j) for j=1,…,2N, we choose one which minimizes the sum of the absolute values of the left-hand sides of (26) for m=N,…,2N-1 given by
(35)∑m=N2N-1‍|σTHN,mσ|.
In other words, to obtain the theoretically unique solution of (26) for m=0,…,2N-1, first we find 2N candidates from (26) for m=0,…,N-1 and then choose the one which best explains the remaining equations (26) for m=N,…,2N-1.


### 3.2. Reconstruction ofμk, νk, and N

To estimateN, following the method for the dipole source model [10], we assume that there are N′(>N) dipole-quadrupole sources and then estimate ζk for k=1,2,…,N′ using the algorithm in Section 3.1. Once they are obtained, μk and νk for k=1,2,…,N′ are linearly solved using (24) for n=0,1,…,2N′-1. Then, we compute |μn+1/μn| and |νn+1/νn| for n=1,2,…,N′-1, which are expected to be sufficiently small when n=N. Practically, we estimate N 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 thez-Coordinates

To determine thez-coordinates of the dipole-quadrupole positions, after we determine the xy-coordinates of the source, we use the method proposed in [7]. We put l=n+1 and m=n in (16) and (17), where n∈ℕ*. From the identity rn+1Y^n+1,n=(2n+1)!(x+iy)nz [23], we have
(36)∑k=1N‍((∇(∇(x+iy)nz×r))T|r=rk∇(x+iy)nz|r=rk·(rk×pk)∑k=1N‍+Qk:(∇(∇(x+iy)nz×r))T|r=rk)=2n+3μ0∫Γ‍n·B(r)(x+iy)nzdS.
Now, let
(37)ϕj(t)≡(t-ζj)2∏l=1,l≠jN‍(t-ζl)3≡∑n=03N-1‍αntn.
Then, from ∑n=03N-1‍αn × (18) and ∑n=03N-1‍αn × (36), we obtain
(38)∑k=1N‍((∇(∇ϕj(x+iy)×r))T|r=rk∇ϕj(x+iy)|r=rk·(rk×pk)+Qk:(∇(∇ϕj(x+iy)×r))T|r=rk)=1μ0∑n=03N-1‍(2n+1)αn∫Γ‍n·B(r)(x+iy)ndS,(39)∑k=1N‍((∇(∇ϕj(x+iy)z×r))T|r=rk∇ϕj(x+iy)z|r=rk·(rk×pk)+Qk:(∇(∇ϕj(x+iy)z×r))T|r=rk)=1μ0∑n=03N-1‍(2n+3)αn∫Γ‍n·B(r)(x+iy)nzdS,
respectively. In the left-hand sides of (38) and (39),
(40)∇ϕj(x+iy)|r=rk=(ϕj′(ζk)iϕj′(ζk)0)=0,(41)∇ϕj(x+iy)z|r=rk=(zϕj′(ζk)izϕj′(ζk)ϕj(ζk))=0
since from (37) it holds that ϕj(ζk)=ϕj′(ζk)=0 for all k=1,2,…,N irrespective of whether or not k is equal to j. Also, we have
(42)∇(∇ϕj(x+iy)×r)|r=rk=δjkΦj,(43)∇(∇ϕj(x+iy)z×r)|r=rk=zjδjkΦj,
where δjk is the Kronecker delta, and
(44)Φj=2∏l=1,l≠jN‍(ζj-ζl)3(izj-zj-iζj-zj-izjζj000).
See Appendix D for the derivation of (42) and (43).Substituting (40) through (43) into (38) and (39) gives
(45)Qj:ΦjT=1μ0∑n=03N-1‍(2n+1)αn∫Γ‍n·B(r)(x+iy)ndS,Qj:zjΦjT=1μ0∑n=03N-1‍(2n+3)αn∫Γ‍n·B(r)(x+iy)nzdS,
from which we obtain
(46)zj=∑n=03N-1(2n+3)αn∫Γ‍n·B(r)(x+iy)nzdS∑n=03N-1(2n+1)αn∫Γ‍n·B(r)(x+iy)ndS.
This is the reconstruction formula of the z-coordinates: once ζk for k=1,2,…,N are obtained, then zj are determined using (46) for each j=1,2,…,N with αn 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 areN′(>N) dipole-quadrupole sources. Using the algorithm in Step 2 (1)–(5) below, where N is replaced by N′, obtain μk and νk for k=1,2,…,N′.
(2) Compute|μn+1/μn| and |νn+1/νn| for n=1,2,…,N′-1. When |μN^+1/μN^|<ϵμ and |νN^+1/νN^|<ϵν for the thresholds ϵμ and ϵν, which are appropriately set a priori, we estimate N to be N^.Step 2 (estimate the position parameters ofN dipole-quadrupole sources).
(1) Givenn·B on Γ, compute cn in (23) for n=0,1,…,4N-1 and construct the Hankel matrices (27) for m=0,1,…,2N-1.
(2) Solve the simultaneous second-degree equation (26) for m=0,1,…,N-1, using the Gröbner bases to get 2N solutions.
(3) Out of the2N solutions, choose one solution σ1,…,σN that minimizes (35).
(4) SolvezN+σ1zN-1+⋯+σN=0 to obtain ζ1,…,ζN.
(5) Solve (24) for μk and νk, where k=1,…,N.
(6) Using the obtainedζ1,…,ζN, compute αn in (37) for n=0,1,…,3N-1 and determine z1,…,zN using (46).


## 3.1. Algebraic Reconstruction of thex- and y-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 l=m=n, where n∈ℕ. From the relationship rnYnn(θ,ϕ)=(2n-1)!!(x+iy)n [23], we have
(18)∑k=1N‍((∇(∇(x+iy)n×r))T|r=rk∇(x+iy)n|r=rk·(rk×pk)∑k=1N‍+Qk:(∇(∇(x+iy)n×r))T|r=rk)=2n+1μ0∫Γ‍n·B(r)(x+iy)ndS.
The primes in both sides are removed for simplicity. Since it holds that(19)∇(x+iy)n=n(1i0)(x+iy)n-1,(∇(∇(x+iy)n×r))T=n(i(n-1)z(x+iy)n-2-(n-1)z(x+iy)n-2i(x+iy)n-1-(n-1)z(x+iy)n-2-i(n-1)z(x+iy)n-2-(x+iy)n-1-in(x+iy)n-1n(x+iy)n-10),we obtain
(20)∇(x+iy)|r=rk·(rk×pk)=n[rk×pk]x+iy(xk+iyk)n-1,(21)Qk:(∇(∇(x+iy)n×r))T|r=rk=n(n-1)((i(qxx,k-qyy,k)-2qxy,k)zk=n(n-1)+(qyz,k-iqxz,k)(xk+iyk))(xk+iyk)n-2,
where [rk×pk]x+iy represents (x-component of rk×pk) + i×(y-component of rk×pk). See Appendix C for the derivation of (21). Substituting (20) and (21) into (18) gives
(22)∑k=1N‍(+(qyz,k-iqxz,k)(xk+iyk))(xk+iyk)n-2[rk×pk]x+iy(xk+iyk)n-1+(n-1)×((i(qxx,k-qyy,k)-2qxy,k)zk+(qyz,k-iqxz,k)(xk+iyk))(xk+iyk)n-2)=2n+1μ0n∫Γ‍n·B(x+iy)ndS.
We now put
(23)ζk=xk+iyk,μk=[rk×pk]x+iy,νk=(i(qxx,k-qyy,k)-2qxy,k)zk+(qyz,k-iqxz,k)(xk+iyk),cn=2n+3μ0(n+1)∫Γ‍n·B(x+iy)n+1dS
and rewrite (22) as
(24)∑k=1N‍μkζkn+n∑k=1N‍νkζkn-1=cn(n∈ℕ*).
Comparing with the special case where the dipoles were distributed in a plane parallel to the xy-plane in [20], one finds that νk has an extra term: (qyz,k-iqxz,k)(xk+iyk).Equation (24) has the same form as that of (6) in [17], and consequently the dipole-quadrupole position projected on the xy-plane, ζk(k=1,2,…,N), can be algebraically reconstructed from cn(n=0,1,…,2N-1), which is shown as follows. In Section 3.1, we assume that N is known a priori and construct an algorithm to estimate ζk in (24). A method to estimate N as well as μk and νk is then described in Section 3.2.Let us defineσ1,σ2,…,σN by the coefficients of
(25)∏k=1N‍(ζ-ζk)=ζN+σ1ζN-1+⋯+σN.
It was proved in [17] that (24) can be transformed into the second-degree equations for σ1,σ2,…,σN:
(26)σTHN,mσ=0(m∈ℕ*),
where σ=(σN,…,σ1,1)T and HN,m is the (N+1)×(N+1) Hankel matrix given by
(27)HN,m=(cmcm+1⋯cm+Ncm+1cm+2⋯cm+N+1⋮⋮⋮⋮cm+Ncm+N+1⋯cm+2N).
Moreover, it has been shown that the 2N equations (26) for m=0,1,…,2N-1 can be turned into linear equations as follows. First, from those 2N equations, we make 2N-1 equations
(28)σTHN,m(1)σ=0(m=0,1,…,2N-2),
where
(29)HN,m(1)≡cm+1HN,m-cmHN,m+1.
By this definition, HN,m(1) are the Hankel matrices whose (1,1)-components are zero. We denote their components by
(30)HN,m(1)=(0cm+1(1)⋯cm+N(1)cm+1(1)cm+2(1)⋯cm+N+1(1)⋮⋮⋮⋮cm+N(1)cm+N+1(1)⋯cm+2N(1)).
Next, from those 2N-1 equations, we make 2N-2 equations:
(31)σTHN,m(2)σ=0(m=0,1,…,2N-3),
where
(32)HN,m(2)=cm+2(1)HN,m(1)-cm+1(1)HN,m+1(1).
By definition, HN,m(2) are the Hankel matrices expressed by
(33)HN,m(2)=(00cm+2(2)⋯cm+N(2)0cm+2(2)⋯⋯cm+N+1(2)cm+2(2)⋯⋯⋯cm+N+1(2)⋮⋮⋮⋮⋮cm+N(2)cm+N+1(2)⋯⋯cm+2N(2)).
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 jth line in the anti-diagonal direction for j=3,…,2N-1, we obtain
(34)σTHN,m(j)σ=0(m=0,1,…,2N-1-j),
where HN,m(j) are the Hankel matrices whose components along the first through jth lines in the anti-diagonal direction are zero. Here, (34) for j=2N-1 is a linear equation for σN. Solving it and substituting the solution into (34) for j=2N-2 gives a linear equation for σN-1. Repeating the backward substitution gives σN-2,…,σ1, successively. The conditions for the unique solvability of these linear equations are ζi≠ζj for i≠j and νk≠0. See [17] for more details. Once σ1,…,σN are obtained, we have ζk as the roots of ζN+σ1ζN-1+⋯+σN=0.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 σN,…,σ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 m=0,…,N-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 2N solutions, (σ1,j,…,σN,j) for j=1,…,2N, we choose one which minimizes the sum of the absolute values of the left-hand sides of (26) for m=N,…,2N-1 given by
(35)∑m=N2N-1‍|σTHN,mσ|.
In other words, to obtain the theoretically unique solution of (26) for m=0,…,2N-1, first we find 2N candidates from (26) for m=0,…,N-1 and then choose the one which best explains the remaining equations (26) for m=N,…,2N-1.


## 3.2. Reconstruction ofμk, νk, and N

To estimateN, following the method for the dipole source model [10], we assume that there are N′(>N) dipole-quadrupole sources and then estimate ζk for k=1,2,…,N′ using the algorithm in Section 3.1. Once they are obtained, μk and νk for k=1,2,…,N′ are linearly solved using (24) for n=0,1,…,2N′-1. Then, we compute |μn+1/μn| and |νn+1/νn| for n=1,2,…,N′-1, which are expected to be sufficiently small when n=N. Practically, we estimate N 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 thez-Coordinates

To determine thez-coordinates of the dipole-quadrupole positions, after we determine the xy-coordinates of the source, we use the method proposed in [7]. We put l=n+1 and m=n in (16) and (17), where n∈ℕ*. From the identity rn+1Y^n+1,n=(2n+1)!(x+iy)nz [23], we have
(36)∑k=1N‍((∇(∇(x+iy)nz×r))T|r=rk∇(x+iy)nz|r=rk·(rk×pk)∑k=1N‍+Qk:(∇(∇(x+iy)nz×r))T|r=rk)=2n+3μ0∫Γ‍n·B(r)(x+iy)nzdS.
Now, let
(37)ϕj(t)≡(t-ζj)2∏l=1,l≠jN‍(t-ζl)3≡∑n=03N-1‍αntn.
Then, from ∑n=03N-1‍αn × (18) and ∑n=03N-1‍αn × (36), we obtain
(38)∑k=1N‍((∇(∇ϕj(x+iy)×r))T|r=rk∇ϕj(x+iy)|r=rk·(rk×pk)+Qk:(∇(∇ϕj(x+iy)×r))T|r=rk)=1μ0∑n=03N-1‍(2n+1)αn∫Γ‍n·B(r)(x+iy)ndS,(39)∑k=1N‍((∇(∇ϕj(x+iy)z×r))T|r=rk∇ϕj(x+iy)z|r=rk·(rk×pk)+Qk:(∇(∇ϕj(x+iy)z×r))T|r=rk)=1μ0∑n=03N-1‍(2n+3)αn∫Γ‍n·B(r)(x+iy)nzdS,
respectively. In the left-hand sides of (38) and (39),
(40)∇ϕj(x+iy)|r=rk=(ϕj′(ζk)iϕj′(ζk)0)=0,(41)∇ϕj(x+iy)z|r=rk=(zϕj′(ζk)izϕj′(ζk)ϕj(ζk))=0
since from (37) it holds that ϕj(ζk)=ϕj′(ζk)=0 for all k=1,2,…,N irrespective of whether or not k is equal to j. Also, we have
(42)∇(∇ϕj(x+iy)×r)|r=rk=δjkΦj,(43)∇(∇ϕj(x+iy)z×r)|r=rk=zjδjkΦj,
where δjk is the Kronecker delta, and
(44)Φj=2∏l=1,l≠jN‍(ζj-ζl)3(izj-zj-iζj-zj-izjζj000).
See Appendix D for the derivation of (42) and (43).Substituting (40) through (43) into (38) and (39) gives
(45)Qj:ΦjT=1μ0∑n=03N-1‍(2n+1)αn∫Γ‍n·B(r)(x+iy)ndS,Qj:zjΦjT=1μ0∑n=03N-1‍(2n+3)αn∫Γ‍n·B(r)(x+iy)nzdS,
from which we obtain
(46)zj=∑n=03N-1(2n+3)αn∫Γ‍n·B(r)(x+iy)nzdS∑n=03N-1(2n+1)αn∫Γ‍n·B(r)(x+iy)ndS.
This is the reconstruction formula of the z-coordinates: once ζk for k=1,2,…,N are obtained, then zj are determined using (46) for each j=1,2,…,N with αn 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 areN′(>N) dipole-quadrupole sources. Using the algorithm in Step 2 (1)–(5) below, where N is replaced by N′, obtain μk and νk for k=1,2,…,N′.
(2) Compute|μn+1/μn| and |νn+1/νn| for n=1,2,…,N′-1. When |μN^+1/μN^|<ϵμ and |νN^+1/νN^|<ϵν for the thresholds ϵμ and ϵν, which are appropriately set a priori, we estimate N to be N^.Step 2 (estimate the position parameters ofN dipole-quadrupole sources).
(1) Givenn·B on Γ, compute cn in (23) for n=0,1,…,4N-1 and construct the Hankel matrices (27) for m=0,1,…,2N-1.
(2) Solve the simultaneous second-degree equation (26) for m=0,1,…,N-1, using the Gröbner bases to get 2N solutions.
(3) Out of the2N solutions, choose one solution σ1,…,σN that minimizes (35).
(4) SolvezN+σ1zN-1+⋯+σN=0 to obtain ζ1,…,ζN.
(5) Solve (24) for μk and νk, where k=1,…,N.
(6) Using the obtainedζ1,…,ζN, compute αn in (37) for n=0,1,…,3N-1 and determine z1,…,zN using (46).


## 4. Numerical Simulations

In this section, our algorithm is verified numerically.Γ was assumed to be a sphere with R=0.12m 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 areN=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.Table 1
Dipole-quadrupole parameters whenN=2.
k(pole no.)
r
k
p
k
(
q
x
x
,
k
,
q
x
y
,
k
,
q
x
z
,
k
,
q
y
y
,
k
,
q
y
z
,
k
,
q
z
z
,
k
)
1
(0.02, 0.04, 0.07) m
(10, 0, 10) nAm
(1, 1, 0, 2, 1, 0) × 10−1 nAm2
2
(0.02, −0.04, 0.07) m
(10, 0, −10) nAm
(1, 1, 0, −2, 1, 0) × 10−1 nAm2First, to determine the number of the sourcesN, we assume that there are N′=3 dipole-quadrupole sources. Figure 2 shows the obtained |μk| and |νk| for k=1,2,3. One observes that |μ3| and |ν3| are much smaller than |μk| and |νk| for k=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 N=2 poles.Estimation ofN: the obtained |μk| (a) and |νk| (b) for k=1,2,3 when assuming N′=3. Since |μ3| and |ν3| are much smaller than |μk| and |νk| for k=1,2, respectively (|μ3/μ2|=3.8×10-3 and |ν3/ν2|=7.2×10-3), we can judge that there are substantially N=2 dipole-quadrupole sources.
(a)
(b)Figure3 shows the localization result projected on the xy- and yz-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.Localization result whenN=2 where the noise level σ=10-2.
(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 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.Figure 4
Standard deviation of the localization errors with respect to the noise levelσ.


### 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(0.02,±d,0.07) where 2d 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.Figure 5
Relative localization error with respect to the distance between the two dipole-quadrupole sources.


### 4.4. Example WhenN=3

We examine the case when there areN=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 N is known a priori. Figure 6 shows the localization result. Comparing with the result when N=2 in Figure 3, the biases as well as the standard deviations of the localization errors become larger when N=3. It is numerically shown in [26] that, when using dipole-only source model, stability gets worse when N increases and it is difficult to estimate the dipoles greater than about N=5. The example in this subsection shows that it is more difficult to estimate the large number (N≥3) of dipole-quadrupole sources. This is considered due to fast spatial decay of the quadrupolar field. However, the dipole-quadrupole source even with N=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.Table 2
dipole-quadrupole parameters whenN=3.
k(pole no.)
r
k
p
k
(
q
x
x
,
k
,
q
x
y
,
k
,
q
x
z
,
k
,
q
y
y
,
k
,
q
y
z
,
k
,
q
z
z
,
k
)
1
(
0.02
,
-
0.04,0.07
) mm
(
10,0
,
10
) nAm
(
1,1
,
0,2
,
1,0
)
×
1
0
-
1 nAm2
2
(
-
0.04,0
,
0.07
) mm
(
10,0
,
-
10
) nAm
(
1,1
,
0
,
-
2,1
,
0
)
×
1
0
-
1 nAm2
3
(
0,0.05,0.07
) mm
(
0,10,10
) nAm
(
2,1
,
0,1
,
1,0
)
×
1
0
-
1 nAm2Localization result whenN=3.
(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 ofr=5mm and a height of h=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.Figure 7
Source of distributed dipoles on a half-cylinder.Case 1.
A single half-cylinder source atr1=(20,50,20) mm. The vectors which determine the posture of the cylinder, e12 and e11, are set to be er and eθ, respectively, where er and eθ are the unit vectors in the r and θ directions at r1. (see Figure 8(a)).  In this case, the total dipole moment p1 is almost parallel to r1; 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.Cases1 and 2 of the cylindrically distributed dipoles.
(a)
(b)Case 2.
The same half-cylinder atr1 as in Case 1, but with e12=-eθ and e11=er (Figure 8(b)). In this case, the angle between r1 and p1 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. N′=3 was assumed to estimate N.


#### 4.5.1. Case 1

Figure9 shows the estimated |μk| normalized by max{μk} using DM. From the fact that |μ3/μ2|=0.005, N 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 |μk| normalized by max{μk} and |νk| normalized by max{νk}, respectively, using DQM. From |μ3/μ2|=0.041 and |ν3/ν2|=0.004, N is estimated to be one.Figure 9
Estimation ofN in Case 1 using the dipole model (DM) when assuming N′=3. From |μ3/μ2|=0.005, DM estimates N=2.Estimation ofN in Case 1 using the dipole-quadrupole model (DQM) when assuming N′=3. From |μ3/μ2|=0.041 and |ν3/ν2|=0.004, DQM estimates N=1.
(a)
(b)Figure11 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.Localization result of the cylindrically distributed dipole sources using dipole model (DM) and dipole-quadrupole model (DQM).
(a)
(b)


#### 4.5.2. Case 2

Figures12, 13(a), and 13(b) show the estimated |μk| using DM and |μk| and |νk| using DQM, respectively. In this Case, from |μ2/μ1|=0.006 in Figure 12 using DM and |μ2/μ1|=0.006 and |ν2/ν1|=0.073 in Figure 13 using DQM, both model estimate N=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.Figure 12
Estimation ofN in Case 1 using the dipole model (DM) when assuming N′=3. From |μ3/μ2|=0.006, DM estimates N=1.Estimation ofN in Case 1 using the dipole-quadrupole model (DQM) when assuming N′=3. From |μ3/μ2|=0.006 and |ν3/ν2|=0.072, DQM estimates N=1.
(a)
(b)Localization result of the cylindrically distributed dipole sources using dipole model (DM) and dipole-quadrupole model (DQM).
(a)
(b)


## 4.1. Demonstration of Our Method

In this subsection, We examine the case where there areN=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.Table 1
Dipole-quadrupole parameters whenN=2.
k(pole no.)
r
k
p
k
(
q
x
x
,
k
,
q
x
y
,
k
,
q
x
z
,
k
,
q
y
y
,
k
,
q
y
z
,
k
,
q
z
z
,
k
)
1
(0.02, 0.04, 0.07) m
(10, 0, 10) nAm
(1, 1, 0, 2, 1, 0) × 10−1 nAm2
2
(0.02, −0.04, 0.07) m
(10, 0, −10) nAm
(1, 1, 0, −2, 1, 0) × 10−1 nAm2First, to determine the number of the sourcesN, we assume that there are N′=3 dipole-quadrupole sources. Figure 2 shows the obtained |μk| and |νk| for k=1,2,3. One observes that |μ3| and |ν3| are much smaller than |μk| and |νk| for k=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 N=2 poles.Estimation ofN: the obtained |μk| (a) and |νk| (b) for k=1,2,3 when assuming N′=3. Since |μ3| and |ν3| are much smaller than |μk| and |νk| for k=1,2, respectively (|μ3/μ2|=3.8×10-3 and |ν3/ν2|=7.2×10-3), we can judge that there are substantially N=2 dipole-quadrupole sources.
(a)
(b)Figure3 shows the localization result projected on the xy- and yz-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.Localization result whenN=2 where the noise level σ=10-2.
(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 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.Figure 4
Standard deviation of the localization errors with respect to the noise levelσ.


## 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(0.02,±d,0.07) where 2d 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.Figure 5
Relative localization error with respect to the distance between the two dipole-quadrupole sources.


## 4.4. Example WhenN=3

We examine the case when there areN=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 N is known a priori. Figure 6 shows the localization result. Comparing with the result when N=2 in Figure 3, the biases as well as the standard deviations of the localization errors become larger when N=3. It is numerically shown in [26] that, when using dipole-only source model, stability gets worse when N increases and it is difficult to estimate the dipoles greater than about N=5. The example in this subsection shows that it is more difficult to estimate the large number (N≥3) of dipole-quadrupole sources. This is considered due to fast spatial decay of the quadrupolar field. However, the dipole-quadrupole source even with N=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.Table 2
dipole-quadrupole parameters whenN=3.
k(pole no.)
r
k
p
k
(
q
x
x
,
k
,
q
x
y
,
k
,
q
x
z
,
k
,
q
y
y
,
k
,
q
y
z
,
k
,
q
z
z
,
k
)
1
(
0.02
,
-
0.04,0.07
) mm
(
10,0
,
10
) nAm
(
1,1
,
0,2
,
1,0
)
×
1
0
-
1 nAm2
2
(
-
0.04,0
,
0.07
) mm
(
10,0
,
-
10
) nAm
(
1,1
,
0
,
-
2,1
,
0
)
×
1
0
-
1 nAm2
3
(
0,0.05,0.07
) mm
(
0,10,10
) nAm
(
2,1
,
0,1
,
1,0
)
×
1
0
-
1 nAm2Localization result whenN=3.
(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 ofr=5mm and a height of h=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.Figure 7
Source of distributed dipoles on a half-cylinder.Case 1.
A single half-cylinder source atr1=(20,50,20) mm. The vectors which determine the posture of the cylinder, e12 and e11, are set to be er and eθ, respectively, where er and eθ are the unit vectors in the r and θ directions at r1. (see Figure 8(a)).  In this case, the total dipole moment p1 is almost parallel to r1; 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.Cases1 and 2 of the cylindrically distributed dipoles.
(a)
(b)Case 2.
The same half-cylinder atr1 as in Case 1, but with e12=-eθ and e11=er (Figure 8(b)). In this case, the angle between r1 and p1 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. N′=3 was assumed to estimate N.


### 4.5.1. Case 1

Figure9 shows the estimated |μk| normalized by max{μk} using DM. From the fact that |μ3/μ2|=0.005, N 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 |μk| normalized by max{μk} and |νk| normalized by max{νk}, respectively, using DQM. From |μ3/μ2|=0.041 and |ν3/ν2|=0.004, N is estimated to be one.Figure 9
Estimation ofN in Case 1 using the dipole model (DM) when assuming N′=3. From |μ3/μ2|=0.005, DM estimates N=2.Estimation ofN in Case 1 using the dipole-quadrupole model (DQM) when assuming N′=3. From |μ3/μ2|=0.041 and |ν3/ν2|=0.004, DQM estimates N=1.
(a)
(b)Figure11 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.Localization result of the cylindrically distributed dipole sources using dipole model (DM) and dipole-quadrupole model (DQM).
(a)
(b)


### 4.5.2. Case 2

Figures12, 13(a), and 13(b) show the estimated |μk| using DM and |μk| and |νk| using DQM, respectively. In this Case, from |μ2/μ1|=0.006 in Figure 12 using DM and |μ2/μ1|=0.006 and |ν2/ν1|=0.073 in Figure 13 using DQM, both model estimate N=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.Figure 12
Estimation ofN in Case 1 using the dipole model (DM) when assuming N′=3. From |μ3/μ2|=0.006, DM estimates N=1.Estimation ofN in Case 1 using the dipole-quadrupole model (DQM) when assuming N′=3. From |μ3/μ2|=0.006 and |ν3/ν2|=0.072, DQM estimates N=1.
(a)
(b)Localization result of the cylindrically distributed dipole sources using dipole model (DM) and dipole-quadrupole model (DQM).
(a)
(b)


## 4.5.1. Case 1

Figure9 shows the estimated |μk| normalized by max{μk} using DM. From the fact that |μ3/μ2|=0.005, N 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 |μk| normalized by max{μk} and |νk| normalized by max{νk}, respectively, using DQM. From |μ3/μ2|=0.041 and |ν3/ν2|=0.004, N is estimated to be one.Figure 9
Estimation ofN in Case 1 using the dipole model (DM) when assuming N′=3. From |μ3/μ2|=0.005, DM estimates N=2.Estimation ofN in Case 1 using the dipole-quadrupole model (DQM) when assuming N′=3. From |μ3/μ2|=0.041 and |ν3/ν2|=0.004, DQM estimates N=1.
(a)
(b)Figure11 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.Localization result of the cylindrically distributed dipole sources using dipole model (DM) and dipole-quadrupole model (DQM).
(a)
(b)


## 4.5.2. Case 2

Figures12, 13(a), and 13(b) show the estimated |μk| using DM and |μk| and |νk| using DQM, respectively. In this Case, from |μ2/μ1|=0.006 in Figure 12 using DM and |μ2/μ1|=0.006 and |ν2/ν1|=0.073 in Figure 13 using DQM, both model estimate N=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.Figure 12
Estimation ofN in Case 1 using the dipole model (DM) when assuming N′=3. From |μ3/μ2|=0.006, DM estimates N=1.Estimation ofN in Case 1 using the dipole-quadrupole model (DQM) when assuming N′=3. From |μ3/μ2|=0.006 and |ν3/ν2|=0.072, DQM estimates N=1.
(a)
(b)Localization result of the cylindrically distributed dipole sources using dipole model (DM) and dipole-quadrupole model (DQM).
(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 x- and y-coordinates of the dipole-quadrupole source by solving the simultaneous second-degree equations, with the result of which the z-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.



---

*Source: 413980-2013-03-26.xml*