An Explicit Method for Inverse Reconstruction of Equivalent Current Dipoles and Quadrupoles

The purpose of the neuromagnetic inverse problem is to reconstruct primary neural current from measured MEG data. It is known that this inverse problem is ill-posed: uniqueness of the solution to the inverse problem is not guaranteed without a priori assumptions on the current source (Fokas et al. (2004)), and, even when using the source model that can be uniquely reconstructed, the obtained solution changes very sensitively depending on the noise contained in MEG data. Thus, employment of a stable algorithm is highly required in order that MEG becomes a non-invasive brain monitoring tool with not only high temporal resolution but also high spatial resolution. Basically, conventional methods are categorized into two groups: parametric approaches and imaging approaches. See the detailed list of references in Baillet et al. (2001). The former methods assume that the current source can be represented by a relatively small number of equivalent current dipoles (ECDs). This source model is shown to be uniquely reconstructed from radial MEG data, except the radial component of the dipole moment, when the head is assumed to consist of concentric spheres. The usual algorithm for this source model is the non-linear least-squares method that minimizes the squared error of data and the forward solution. An advantage of this algorithm is that the parametric description allows us the accurate estimation of the center position of the activated region. However, the problems of this algorithm are: 1) it requires an initial solution, 2) it requires an iterative computing of forward solution, 3) it is often trapped by the local minimum of the cost function, which gives a solution far from the true one, 4) estimation of the number of ECDs is difficult, and 5) spatial extent of the source is not considered. The secondmethods, imaging approaches, assume that there exist current dipoles at the nodes of artificial meshes on the cerebral surface, and solve a system of linear equations for the dipole moments at the fixed positions. An advantage of this method is that it can describe the spatial distribution of the neural current. However, the problems of this algorithm are: 1) the solution is not unique, 2) adding a regularization term often gives a unique but over-smoothed solution, 3) choice of the regularization parameters strongly affects the solution. Recently, a method with the multipolar representation of the source has been developed that incorporates some of the advantages of the above twomethods, and has attracted considerable attention (Irimia et al. (2009); Jerbi et al. (2002; 2004); Nara (2008a); Nolte et al. (1997; 2000)). In this model, instead of expressing the current source by an equivalent current dipole, 6


Introduction
The purpose of the neuromagnetic inverse problem is to reconstruct primary neural current from measured MEG data. It is known that this inverse problem is ill-posed: uniqueness of the solution to the inverse problem is not guaranteed without ap r i o r iassumptions on the current source (Fokas et al. (2004)), and, even when using the source model that can be uniquely reconstructed, the obtained solution changes very sensitively depending on the noise contained in MEG data. Thus, employment of a stable algorithm is highly required in order that MEG becomes a non-invasive brain monitoring tool with not only high temporal resolution but also high spatial resolution. Basically, conventional methods are categorized into two groups: parametric approaches and imaging approaches. See the detailed list of references in Baillet et al. (2001). The former methods assume that the current source can be represented by a relatively small number of equivalent current dipoles (ECDs). This source model is shown to be uniquely reconstructed from radial MEG data, except the radial component of the dipole moment, when the head is assumed to consist of concentric spheres. The usual algorithm for this source model is the non-linear least-squares method that minimizes the squared error of data and the forward solution. An advantage of this algorithm is that the parametric description allows us the accurate estimation of the center position of the activated region. However, the problems of this algorithm are: 1) it requires an initial solution, 2) it requires an iterative computing of forward solution, 3) it is often trapped by the local minimum of the cost function, which gives a solution far from the true one, 4) estimation of the number of ECDs is difficult, and 5) spatial extent of the source is not considered. The second methods, imaging approaches, assume that there exist current dipoles at the nodes of artificial meshes on the cerebral surface, and solve a system of linear equations for the dipole moments at the fixed positions. An advantage of this method is that it can describe the spatial distribution of the neural current. However, the problems of this algorithm are: 1) the solution is not unique, 2) adding a regularization term often gives a unique but over-smoothed solution, 3) choice of the regularization parameters strongly affects the solution.
Recently, a method with 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 (Irimia et al. (2009) ;Jerbi et al. (2002;; Nara (2008a); Nolte et al. (1997;). In this model, instead of expressing the current source by an equivalent current dipole, 6 www.intechopen.com an equivalent dipole and quadrupole (Jerbi et al. (2002;) or an equivalent dipole and octopole (Nolte et al. (1997;) are used where the quadrupole or octopole is determined depending on the spatial extent of the support of the current source. Jerbi et al. (2004) 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. We considered a two-dimensional (2D) problem using a complex analysis framework and proposed an explicit method to reconstruct the dipole and quadrupole parameters directly from the boundary measurements of the meromorphic function (Nara (2008a)). In Nara (2008b), we proposed an explicit method for 3D case when the dipoles were distributed in a plane parallel to the xy-plane. The aim of this chapter is to describe explicit methods for the equivalent current dipole model and the equivalent current dipole-quadrupole model in Nara et al. (2007), Nara (2008a), andNara (2008b) and compare them using numerical simulations. Here, the term 'explicit' means that we have an analytic and explicit form of the solution to the inverse problem. As a consequence, the explicit method for the ECD model can resolve the problems 1)-4) in the parametric approach with the conventional, non-linear least-squares method listed above. That is, without an initial solution and iterative computation of forward solution, the algorithm can reconstruct the ECD parameters. The number of the ECDs can be also estimated. With the explicit reconstruction formula, the sensitivity analysis can be conducted in which the estimation error is evaluated with the noise level. From a practical viewpoint, the solution obtained by the explicit method can be used as a good initial solution close to the true one for the iterative non-linear squares methods. Moreover, in order to resolve the problem 5) in the parametric approach, we show an explicit method for the dipole-quadrupole source model. The rest of this chapter is organized as follows. In section 2, we introduce the equivalent current dipole and quadrupole source model, and show how it expresses the spatial extent of the current source. In section 3, an explicit method is shown: subsection 3.1 describes an algorithm for the equivalent current dipole model, whereas subsection 3.2 explains a method for the equivalent current dipole-quadrupole source model. In section 4, both the algorithms are compared with numerical simulations.

Problem setting
In this chapter, we explain our explicit algorithm using the spherical head 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. Here, Ω 1 , Ω 2 /Ω 1 ,a n dΩ 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 modelaswellasthesphericalsensorsurface,themethodcanbeextendedtoamorerealistic case when the head is modeled by a piecewise homogeneous layered domain and the sensors are set on an arbitrarily shaped surface (Nara et al. (2007)). Let us assume that the neural current J p is supported on several domains D k (⊂ Ω 1 ) for k = 1, 2, ··· , N. The 'center' position of D k is denoted by r k =( x k , y k , z k ) T ;t h isist h ema in parameter to be reconstructed. We express this source equivalently by the current dipoles and quadrupoles: is the equivalent current dipole for the source in D k ,and q xx,k q xy,k q xz,k q yx,k q yy,k q yz,k q zx,k q zy,k q zz,k ⎞ ⎠ is the equivalent current quadrupole for the source in D k .N o t eh e r et h a tp k does not depend on the size of D k , while Q k depends on the spatial extent of J p around r k in D k . Q k is a 3 × 3 tensor of order 2 and is called the quadrupole moment tensor. In this chapter, we assume that Q k are of the form ⎛ ⎝ q xx,k q xy,k 0 q yx,k q yy,k 0 00 0 ⎞ ⎠ .
In other words, the quadrupoles are parallel to the xy-plane. Extention to general case where all the components of Q k are not zero is an important aspect of further studies. The validity of the source model (1) is confirmed as follows. Substituting Eq. (1) into the expression of radial MEG given by the Biot-Savart law where μ 0 is the magnetic permeability assumed to be constant in the whole space, we have Here, ':' represents the tensor contraction defined by A : B = ∑ 3 i=1 ∑ 3 j=1 A ij B ij for second order tensors A and B whose (i, j) components are given by A ij and B ij , respectively, and X r is the cross product tensor defined by and hence is written as  Jerbi et al. (2002) (when N = 1), which was obtained by truncating the Taylor series expansion of radial MEG, up to the secondly dominant terms, generated by J p with spatial extent. The first term in Eq. (3) is the magnetic field created by the equivalent current dipole p k , and the second term is the magnetic field created by the equivalent current quadrupole Q k . Hence we call Eq. (1) the equivalent current dipole-quadrupole source model. Using this source model (1), our inverse problem is formulated as follows: reconstruct r k , p k , Q k and N from radial MEG data on Γ. Since estimation of p k and Q k is a linear problem if r k and N are determined, we restrict our interest in this chapter to estimation of r k and N.

Explicit method
In this section, we show an explicit method for the source model (1). In subsection 3.1, a method for the dipole source model (when Q k = 0) is shown. An algorithm for the dipole-quadrupole source model is described in subsection 3.2. Both the algorithms are based on the multipole expansion of radial MEG: Eq. (2) can be expressed by the multipole expansion and P m l (cos θ) are the associated Legendre polynomials. As shown in Eq. (84) in Jerbi et al. (2002), the multipole moment has the following relationship with J p : Substituting Eq. (1) into Eq. (4) and using Eq. (17) 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 Γ (Alvarez (1991)): where n ′ is the outward unit normal vector at r ′ on Γ. Equating Eqs. (5) and (6) for l ≥ m ≥ 0 gives algebraic equations relating the unknown parameters of the sources to the radial MEG data.

Explicit method for the equivalent current dipoles
First, we describe an explicit algorithm for the equivalent current dipole source model (Q k = 0). Let us consider Eqs. (7) for l = m, which are called 'sectorial harmonics' components. Since in this case it holds that (Hobson (1931)) The prime characters in the integrands have been omitted for brevity. Define now the following quantities: Eq. (8) can be rewritten as where The problem of determining S k , μ k and N in Eqs. (10) from c m (m = 0, 1, ···) appears in many inverse problems such as computed tomography (Golub et al. (1997)), EEG inversion (El-Badia et al. (2000); Nara et al. (2003)), MEG inversion (Nara et al. (2007)), and locating the zeros of analytic functions (Kravanja et al. (1994)). To this problem, an algebraic algorithm called Prony's method and its modified/extended algorithms have been proposed that enable us to reconstruct the source parameter S k , μ k and N from c m . The essence of Prony's method is as follows (See e.g. Elad et al. (2004)). First, let us assume that N is fixed and known apriori. Estimation of N is explained afterwards. For N positions S 1 , ··· , S N ,letusdefineσ 1 , ··· , σ N by the coefficients of the polynomial Then, the following difference equations holds Thus, from the 2N linear equations (11) Here, since the coefficient matrix in Eq. (12) is diagonalized as Hence, when μ k = 0andS i = S j , Eqs. (12) can be uniquely solved for σ 1 , ··· , σ N . Then, S n are obtained as solutions to the Nth degree equation Although it is known that Prony's method is numerically unstable, modified versions of Prony's method have been proposed by using the truncated singular value decomposition (Kumaresan et al. (1982)), by the projection method (Cadzow (1988)), or by repeating Prony's method while changing the model order N (Barone et al. (1998)). Methods to transform Eqs. (10) to an eigenvalue problem have also been proposed (Hua et al. (1990), Luk et al. (1997), Golub et al. (1997), El-Badia et al. (2000)). See Elad et al. (2004) for a review. Onece S k are determined, μ k can be obtained by solving the linear equations (10).
To determine z k we use Eqs. (7) for l = m + 1. Since it holds that (Hobson (1931)) where Hence, after S k and μ k are obtained, z k and the z-component of µ k denoted by [µ k ] z can be determined by solving Eqs. (13). In order to estimate N, we assume that there exist N ′ (> N) dipoles and estimate μ 1 , ··· , μ N ′ using the algorithm mentioned above. Then it would be expected, when the data includes noise, that μ k of the 'spurious' sources k = N + 1, ··· , N ′ aremuchsmallerthanthoseofthe true sources, from which N can be estimated.

Explicit method for the equivalent current dipoles and quadrupoles
Now we describe an explicit algorithm for the equivalent current dipole and quadrupole source model (Q k = 0). Considering again Eqs. (7) for l = m,wehave It is shown in Nara (2008b) that Eq. (14) is rewritten as where ν k =( i(q xx,k − q yy,k ) − (q xy,k + q yx,k ))z k .
Eq. (15) has the same form as that of Eq. (6) in Nara (2008a), and can be transformed into the simultaneous second degree equations where Furthermore, the second degree equations (16) for m = 0, 1, ··· ,2N − 1 can be turned into linear equations for σ 1 , ··· , σ N . (See Nara (2008a)). An example of the algorithm when N = 2 is illustrated in Fig. 1: in Step 1), starting with σ T H N,m σ = 0form = 0, 1, 2, 3, we compute their linear combinations to have three equations with the Hankel matrices whose (1,1)-components are zero. In Step 2), the linear combinations of those three equations are computed so that we have two equations with the Hankel matrices whose first and second anti-diagonal components are zero. In Step 3), by computing the linear combinations of those two equations, we arrive at a single equation, which is a linear equation for σ 1 . In Step 4), by substituting the obtained σ 1 into the equations obtained in Step 2), we have linear equations for σ 2 . This is an example how the algorithm proceeds when N = 2. See Nara (2008a) for the detailed algorithm for general N.
Once σ 1 , ··· , σ N are obtained, S 1 , ··· , S N are obtained by solving ζ N + σ 1 ζ N−1 + ···+ σ N = 0. μ k and ν k for k = 1, 2, ··· , N are linearly solved using Eqs. (15) for m = 0, 1, ··· ,2N − 1. Then, we use Eqs. (7) for l = m + 1 which leads to the second degree equations for z k . To estimate N, following the method for the dipole source model, we assume that there are N ′ (> N) dipole-quadrupole sources and then estimate S k as well as μ k and ν k for k = 1, 2, ··· , N ′ .Th enw ec omput e|μ k+1 /μ k | and |ν k+1 /ν k | for k = 1, 2, ··· , N ′ − 1, which are expected to be sufficiently small when k = N. Practically, we estimate N such that these ratios become smaller than some thresholds set apriori. 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 (Barone et al. (1998)  is included in both A and B. When the data includes noise, we choose an element in A and B,s a yS (1) and S (3) , such that the distance between them, |S (1) − S (3) |, is smaller than the distance between the other pair, |S (2) − S (4) |, and estimate the projected position S by S (1) +S (3) 2 . In the simulations in section 4, we use this algorithm.

Numerical simulations
In this section, we compare the explicit method assuming the dipole-quadrupole model (DQM) with the explicit method assuming the dipole model (DM). To model dipoles on cerebral convolutions, we assume that dipoles are placed on a mesh on a half cylinder with ar a d i u so fr = 5mm and a height of h = 5 mm, as shown in Fig. 2. 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 convolutions (Hämäläinen et al. (1993)). We examined the following three cases: • case (i) A single half-cylinder source at r 1 = 70(sin θ 1 cos φ 1 ,sinθ 1 sin φ 1 ,cosθ 1 ) mm, where θ 1 = 10 180 π and φ 1 = π. The vectors which determine the posture of the cylinder, e 12 and e 11 ,aresettobe(0, 0, 1) and (1, 0, 0), respectively. See Fig. 3 left. In this case, the total dipole moment p 1 is nearly parallel to r 1 ; the angle between them is 9.8 degrees. Since a radial dipole is a silent source for radial MEG (Hämäläinen et al. (1993)), this cylindrical source is regarded as being almost quadrupolar.
• case (iii) Two half-cylinder sources of cases (i) and (ii). See Fig. 3 right. The source at r 1 is almost quadrupolar and that at r 2 is a dipole-quadrupole source. Fig. 2. 'Cylindrical source': distributed dipoles on a half cylinder modeling the cerebral convolutions. Fig. 3. Case (i) a single cylindrical source where r 1 is nearly parallel to the equivalent current dipole p 1 (the angle between them is 9.8 degrees); this cylindrical source is almost quadrupolar. Case (ii) a single cylindrical source where r 1 is nearly perpendicular to p 1 (the angle between them is 78 degrees); this cylindrical source has a detectable equivalent dipole moment as well as the equivalent quadrupole moment. Case (iii) two cylindrical sources (combination of cases (i) and (ii)).
We computed the forward solution generated by the 30 (case (i) and (ii)) or 60 (case (iii)) elemental dipoles using Eq. (2). Note that Eq. (3) was not used to compute the theoretical data. The radius of a head was set to be 100 mm. We assumed that the radial component of the magnetic field was measured at M = 361 points distributed uniformly on a sphere with R = 120 mm using the spherical t-design (Saff97 et al. (1997)); it is a set of M points on Γ such that the integral of any polynomial of degree t or less over Γ is equal to the average value of the polynomial over the set of M points. We used M =( t + 1) 2 = 361 points (t = 18) given by Chen et al. (2009). Based on this property, for numerical integration, the integrand values on the N points were summed with equal weights 4π M .
To validate our algorithm for the dipole-quadrupole model, we assumed that the data was available on the whole sphere which enclosed the source. Simulations using sensors on a part of the sphere is left for future studies. To this end, the so-called Signal Space Separation method proposed by Taulu et al. (2005) or the method for stable data continuation from data on the upper hemisphere to those on the lower hemisphere proposed by Popov (2002) could be very useful. Gaussian noise was added to the theoretical forward solution, where the noise level defined by the ratio of the standard deviation of the noise to the root mean squares of the data was 5%. 10 data sets with the different noise added were used for reconstruction. First we show the reconstruction result when N is know ap r i o r i . Table 1 shows the error between the true position and the mean estimated position using 10 data sets. The true and estimated positions projected on the xy-plane are depicted in Fig. 4. We observe that in case (i) the method using DM cannot estimate the source accurately, while the method using DQM can estimate it within an error of 2 mm. This is because the source is almost quadrupolar.
In case (ii), the result using DM is better than that using DQM. In case (iii), the maximum error about the xyp r o j e c t e dp o s i t i o n su s i n gD Q Mi s7 . 6m m , whereas that using DM is 110 mm (for the almost quadrupolar source). However, the z-coordinates is not accurately obtained even when using DQM. This is because the z-coordinates are estimated using the obtained x-andy-coordinates.
Next, we show the case when N is not known apriori. Fig. 5 shows the reconstruction result when assuming N ′ = 2 in cases (i) and (ii) and N ′ = 3 in case (iii). In case (i), when using DQM where N ′ = 2, the two positions are estimated: one is close to the true one, and another is far from the true one (In Fig. 5 top left, the estimated position far from the true position is not seen, since it is out of the figure.) Numbering them 1 and 2, we have |μ 2 /μ 1 | = 5.3e − 4 and |ν 2 /ν 1 | = 9.9e − 4. From this, we can reasonably judge that the second source is spurious due to the noise, and there is a single dipole-quadrupole source. In contrast, when using DM where N ′ = 2, |μ 2 /μ 1 | = 4.3e − 1; μ 1 and μ 2 are almost the same order, and hence we judge that there are two dipoles. The estimated dipoles are close to the side walls of the half cylinder. However, the distance between them is larger than the diameter of the cylinder as in Fig. 5 top left.
In case (ii), when using DQM where N ′ = 2, |μ 2 /μ 1 | = 9.7e − 4and|ν 2 /ν 1 | = 6.5e − 3. Also, when using DM where N ′ = 2, |μ 2 /μ 1 | = 4.1e − 3. Hence, both DQM and DM can estimate the number of the sources. In case (iii), when using DQM where N ′ = 3, |μ 3 /μ 2 | = 9.5e − 4and|ν 3 /ν 2 | = 9.0e − 4, from whichwecanjudgethatN = 2. In contrast, when using DM where N ′ = 3, |μ 3 /μ 2 | = 1.1e − 1; the third source is not much smaller than the second one. Also we observe that, for almost the quadrupole source, the estimated two dipoles are too separated. Fig. 4. The reconstruction result projected on the xy-plane. In this figure, N is assumed to be known apriori. (top left) In case (i), the equivalent dipole is almost directed to the radial direction, and the source is almost quadrupolar. As a result, when using DM with N = 1, the cylindrical source is not at all localized, while DQM with N = 1 well estimates the center of the source. (top right) In case (ii), the equivalent dipole is almost perpendicular to the radial direction. In this case, both DM and DQM can well localize the source. (bottom) In case (iii), both the sources are well estimated when using DQM, while almost the quadrupole source is not at all estimated when using DM. when assuming N ′ = 2. When using DQM, the estimated source far from the true position (which is out of the figure and is not depicted) has much smaller dipole and quadrupole moments than the estimated source close to the true position. In fact, |μ 2 /μ 1 | = 5.3e − 4and |ν 2 /ν 1 | = 9.9e − 4, from which we can judge N = 1. In contrast, when using DM, |μ 2 /μ 1 | = 4.3e − 1 from which we judge that there are two dipoles. Although two dipoles are estimated close to the side walls of the cylindrical surface, the distance between them is larger than the diameter of the cylinder. (top right) In case (ii), |μ 2 /μ 1 | = 9.7e − 4and |ν 2 /ν 1 | = 6.5e − 3 when using DQM where N ′ = 2. Also, |μ 2 /μ 1 | = 4.1e − 3 when using DM. Hence, both DQM and DM can estimate the number of the sources. (bottom) In case (iii), |μ 3 /μ 2 | = 9.5e − 4and|ν 3 /ν 2 | = 9.0e − 4 when using DQM where N ′ = 3, from which we can judge that N = 2. In contrast, when using DM where N ′ = 3, |μ 3 /μ 2 | = 1.1e − 1; the third source is not much smaller than the second one. Also we observe that, for almost the quadrupole source, the estimated two dipoles are too separated.

Conclusion
In this chapter, we introduced the equivalent current dipole-quadrupole source model which has a potential to parametrically represent the spatial extent of the neural current in MEG inverse problem. Then, explicit methods for the equivalent dipole-quadrupole source model as well as the equivalent dipole source model were shown, that enables us to reconstruct the dipole-quadrupole parameters explicitly with MEG data. In numerical simulations, it was suggested that the dipole-quadrupole source model would be useful especially when the elemental dipoles are distributed on the surface of a half cylinder modeling the cerebral convolution such that the equivalent dipole is parallel to the radial direction.

Appendix
It is easy to obtain the dipole terms. For the quadrupole terms, we use the identity for an arbitrary vector field a(r ′ )=(a x (r ′ ), a y (r ′ ), a z (r ′ )) T ,whereT represents the transpose. When inserting the quadrupole terms in Eq. (1) into Eq.
(2), we have from Eq. (17) Here, it holds that and hence (∇ ′ (∇ ′ 1 |r − r ′ | × r ′ )) T | r ′ =r k = 3(r × r k )(r − r k ) |r − r k | 5 + X r |r − r k | 3 . This is a practical book on MEG that covers a wide range of topics. The book begins with a series of reviews on the use of MEG for clinical applications, the study of cognitive functions in various diseases, and one chapter focusing specifically on studies of memory with MEG. There are sections with chapters that describe source localization issues, the use of beamformers and dipole source methods, as well as phase-based analyses, and a step-by-step guide to using dipoles for epilepsy spike analyses. The book ends with a section describing new innovations in MEG systems, namely an on-line real-time MEG data acquisition system, novel applications for MEG research, and a proposal for a helium re-circulation system. With such breadth of topics, there will be a chapter that is of interest to every MEG researcher or clinician.