Brain Topography, Volume 5, Number 2, 1992
95
Magnetoencephalography: A Tool For Functional Brain Imaging Matti S. H~m~Jl~inen* Summary: At present, one of the most promising windowsto the functionalorganizationof the human brain is magnetoencephalography(MEG). Bymappingthe magneticfielddistributionoutsidethe head the sites of neuraleventscan be locatedwith an accuracyof a few millimetersand the temporalevolutionofthe activationcanbe tracedwitha millisecondresolution.Thispaper reviewssomeforwardfieldcalculationapproachessuitable forthe interpretationofthe brain's electromagneticsignals.Inversemodellingwithmultipledipolesis describedin detail.An exampleof the analysis ofthe sornatosensoryevoked-responsesillustratesthepotentialofmultiplesignalclassification(MUSIC)algorithmin findingoptimaldipolepositions. Keywords: Magnetoencephalography;Functionalbrain imaging;Forwardproblem; Currentsource estimation.
Introduction Modern imaging techniques provide high-resolution windows through which the structure and metabolism of the h u m a n brain can be studied in vivo. Computer-assisted X-ray tomography (CT) and magnetic-resonance imaging (MRI) produce high-quality but static maps of the brain's anatomy. Functional information about the brain on a time scale of minutes can be obtained with regional blood flow (RCBF) and positron-emission tomography (PET) measurements, based on changes in cell metabolism as indicated by radioactive tracers. There are only two totally noninvasive windows to study the actual brain function on a millisecond scale: electroencephalography (EEG) and magnetoencephalography (MEG). EEG, employed for about 60 years, involves the measurement of the electric potential distribution by means of electrodes attached on the scalp. In MEG, the minute (50 - 500fT) magnetic fields, produced by the same brain events as in EEG, are detected outside the head with SQUID (Superconducting QUantum Interference Device) magnetometers. The processes which can be studied by means of MEG *LowTemperatureLaboratory,HelsinkiUniversityof Technology, Espoo,Finland. Acceptedfor publication: June23, 1992. I thank Prof. Riitta Hari, Prof. Olli V. Lounasmaa,and Dr. Juha Simola for valuable comments on this manuscript. This work was supported by Finnish National Fund for Researchand Development (SITRA)and InstrumentariumCorp. Correspondenceand reprint requests shouldbe addressedto Matti S. Hamalainen,Low Temperature Laboratory,HelsinkiUniversityof Technology,02150Espoo,Finland. Copyright© 1992Human SciencesPress, Inc.
may be divided into two main categories: evoked responses and spontaneous brain activity. In an evoked-field experiment, sensory stimuli are presented to the subject and the ensuing magnetic field is recorded outside the head at several locations; at each site the stimulus is repeated 50 - 200 times in order to remove incoherent noise by signal averaging. Spontaneous phenomena, both normal and abnormal, resulting from localized events in the brain, e.g., the 10 Hz alpha-rhythm and epileptic activity, can also be investigated provided that the signal-to-noise ratio is sufficient, The first MEG measurements were made with induction coil magnetometers in the late 1960's. At present, multichannel SQUID devices of considerable sophistication are used by several research groups around the world. These instruments include the Helsinki 24-channel system (Kajola et al. 1989; Ahonen et al. 1991), the Siemens 37-channel 'KrenikonTM' (Hoenig et al. 1991), and the 37-channel 'Magnes TM' of Biomagnetic Technologies, Inc. Although the present instruments allow several types of experiments to be carried out, it has become evident that, in order to study complicated cognitive processes and s p o n t a n e o u s activity in h e a l t h y subjects and patients, it is necessary to record the magnetic field over the whole scalp simultaneously. The first pilot evokedresponse measurements with a 122-channel 'Neuromag122' system have already been carried out in our laboratory. Rather than attempting to cover all aspects of MEG recordings and their interpretation, we will focus on the problem of neuromagnetic source modelling, emphasizing the new requirements due to simultaneous full-head
96
H6m~l~inen
proble~ Forward
Maxwetl*s equations (~
Inverse problem
Figure 1. Schematic illustration of the forward and inverse problems.
coverage. We shall elaborate on the different available models for forward field computation and for the source itself. For in-depth treatises on neuromagnetic measurements and their interpretation an interested reader may wish to consult some of the available review articles (Williamson and Kaufman 1981; Sarvas 1987; Romani and Rossini 1988; Ryh~inen et al. 1989; Hari 1991; Williamson et al. 1991; H~im~il~iinenet al. 1993).
Source Modelling The neuromagnetic inverse problem is to estimate the cerebral current sources underlying a measured distribution of the magnetic field (see figure 1). We shall limit our discussion to this task and will not elaborate on the physiological interpretation, although we shall keep in mind that neuromagnetic fields arise in the brain, which imposes constraints on the sources that the measurements might reveal. Contrary to CT and MRI, the inverse problem of MEG does not have a unique solution (Helrnholtz 1853), even if the magnetic field were precisely known everywhere outside the head. Inclusion of EEG information does not help: we must still restrict the viable source current distributions by choosing a suitable source model. In fact, one can find current distributions that are either magnetically silent (MEG vanishes), electrically silent (EEG vanishes), or both. A simple example of a magnetically silent source that produces an electric field is a radial dipole in a spherically symmetric conductor. The opposite is a current loop, which is electrically silent but produces a magnetic field. Because of the nonuniqueness of the inverse problem it is necessary to confine to a restricted class of source currents. We shall here discuss the problems associated with modelling the neural activity as a collection of point sources, current dipoles. An alternative approach of using distributed source m o d e l l i n g was also i n t r o d u c e d at an early stage (Ham~il~iinen and Ilmoniemi 1984; Ilmoniemi et al. 1985) but, due to technical limitations, it has become useful only after the introduction of the multichannel magnetometers (Kullmann et al. 1989; Ribary et al. 1991).
Forward solutions The solution of the inverse problem necessarily involves the forward field and potential calculations. We shall, therefore, next give a synopsis of the results involved in the computation of the electromagnetic field ensuing from neural current sources. We shall also point out some characteristic differences between MEG and EEG, which emerge from the theoretical results. From a known distribution of currents it is, in principle, straightforward to compute the corresponding distribution of the magnetic field and electric potential using Maxwell's equations. To accomplish the calculation it is necessary to know the spatial distribution of the electrical conductivity, the dielectric constant, and the magnetic permeability. However, for MEG and EEG the calculation is simplified because we can omit the time dependent terms from Maxwell's equations: the quasistatic approximation is valid. Furthermore, we can assume that the magnetic permeabilities of the media in question are equal to that of free space. The remaining problem of field computation is illustrated in figure 1. We start from a known neural current distribution, the primary current J--~,which gives rise to a passive ohmic current, commonly called the volume current, P , extending throughout the head. The volume current is directly related to the electric t~otential gradient by the Ohm s law: jv ( r-~ = - o( ~ vV( r--J,where is the conductivity and V the electric potential; EEG involves the measurement of V on the surface of the head. It is w o r t h n o t i n g that the c o n d u c t i v i t y m a y be anisotropic. Computation of the electric potential is then numerically more complex than in the isotropic case as will be discussed later in this article. Once we know both J--~,and J-~, the magnetic field can be calculated from the Biot-Savart's law:
(1)
MEG: A Tool for Functional Brain Imaging
97
Measurernsnt I
Figure 2. The choice of the origin in the sphere model should be based on the local curvature of the skull's inner surface. The two measurement areas shown require different origins, O1 and 02. The corresponding radii of curvature of the two areas are denoted by R1 and R2, respectively. Typical sources that we would expect to locate with _~easur~ments over the areas shown are denoted by Q1 and Q2. (right) A piecewise homogneous conductor. Each region Gi has uniform conductivity ~i. Unit vectors normal to the surfaces are denoted by nij.
The Sphere Model If the conductor has__a simple symmetrical shape, the computation the V and B can become remarkably simple: we may obtain quickly converging series representations or even an analytic expression. The sphere model is often used in MEG and EEG forward field calculations: the head is assumed to consist of several spherically symmetric layers. It turns out that the magnetic field remains the same, irrespective of the radii and conductivities of the shells. As discussed by de Munck (1989), this is the case even in an anisotropic conductor if there is a difference between the radial and the two tangential conducfivities only. One can derive an analytical formulae for B allowing fast evaluation. The sphere model can be readily applied once we have selected a suitable center of symmetry. Since the skull is a poor conductor, major contribution to the magnetic field comes from the currents flowing within the cranial space. The sphere model should, therefore, be fitted to the inner curvature of the skull, as shown in figure 2. This task is straightforward if MR/images from the subject are available. Otherwise, one must rely on an approximation from the outer shape of the head. From the sphere model one can predict that radial sources are not seen in MEG and that deep sources are strongly attenuated. This means that MEG selectively detects the activity of fissural cortex. Fortunately, many interesting areas, such as the primary auditory, motor, somatosensory, and visual cortices, are located in fissures. Sources in them are~ therefore, easily detectable by MEG. Furthermore, the activity of the convexial cortex,
visible in EEG, is filtered out, whereby the signal-to-noise ratio is improved. The electric potential on the surface of a spherically symmetric conductor is affected both by the conductivity profile and anisotropies in conductivity. For example, if a poorly conducting layer (skull) is placed between two relatively good conductors (brain and scalp), the potential pattern caused by a current source inside will be more widespread than in a homogeneous sphere. Another difference is that both radial and deep sources contribute significantly to V. The computation of V requires the evaluation of a series expansion in spherical harmonics (Arthur and Geselowitz 1970; de Munck 1988), which is more time consuming than the computation of B from simple analytical expressions.
Modelling Conductors of Arbitrary Shape Other simple conductor shapes that can be handled analytically are spheroids (de Munck 1988; Cuffin and Cohen 1977) and eccentric spheres (Meijs and Peters 1987). However, if the shape of the conductor is not restricted we must resort to more general numerical methods for computing V and B. In the quasistatic approximation, the divergence of the total current density ~r--) = J-~ (r--) + j g ( r ' ) vanishes: V. J-~ r--) : O, which leads directly to the Poisson's equation for the electric potential: v. (
vv) = v. pp
(2)
With cy and J~ known, we can solve this differential equation numerically by means of the finite element method (FEM). This approach is very general: we can have any conductivity distribution with anisotropies. We shall focus on another possibility, the boundaryelement method. The head is then assumed to consist of several regions with uniform and isotropic conductivities, see figure 2. Using Green's identity it is now possible to convert the differential equation for V into an integral equation involving V on the boundaries of regions of different conductivity only. For a single homogeneous, arbitrarily shaped conductor with conductivity cy,surrounded by an insulator, the equation reads (Geselowitz 1967):
V ( ~ = 2Vo(r-} + G
(3)
where r' is a point on the bounding surface S, V0( ~ is the potential if the conductor had an infinite extent, and dr2 r-( r-+) is the solid angle subtended at r---boya surface element at r-9'.
98
H6m~l~inen
For a solution of this equation, S is approximated by a net of triangles (Oostendorp and van Oosterom 1989; H~im~ilainen and Sarvas 1989). We can then either assume that V is constant within one triangle or use higher-order approximations to yield more accurate results with smaller computational burden (Urankar 1990; de Munck 1992). The discretization results in a set of linear algebraic equations which can be solved for the unknown potential. Once V is known on S, we can use the Geselowitz' formula, derived from Eq. (1):
N0
v ( ~ ' ) v"
1
tr~_r-+,t
x d ~ r--~')
(4) to compute the magnetic field. In one series of simulations (H~im~il~iinenand Sarvas 1989) comparing the sphere model with more accurate, realistically shaped models, it was found that, except for frontotemporal sources, the sphere provides a sufficient accuracy for practical purposes. When the sphere model is imprecise, a brain-shaped homogeneous model, which assumes that the skull is a perfect insulator, is an excellent approximation. It must be noted that if one likes to predict the EEG potential distribution with the same accuracy as that provided by the homogeneous model in MEG, one must apply a multilayer realistically shaped head model with known conductivities and shapes for the brain, the skull, and the scalp. Much less information is thus needed for the analysis of MEG data.
Multidipole Models Formulation of the Problem
An obvious generalization of the popular singledipole model is an assumption of multiple spatially separated dipolar sources. If the distance between the individual dipoles is sufficiently large ( >3 cm) and their orientations are favorable, the field patterns may show only minor overlap and they can be fitted individually using the single-dipole model. An example of this approach is the separation of activities from the first and s e c o n d s o m a t o s e n s o r y cortices (Hari et al. 1984; Kaukoranta et al. 1986; Hari 1990). Similarly, if the temporal behavior of the dipoles differ, it is often possible to recognize each source separately. However, when the sources overlap both temporally and spatially, we must resort to the multidipole approach to obtain correct results. An effective approach to this modelling problem is to take into account the spatiotemporal course of the signals as a whole, instead of considering each time sample separately. This method was first
applied to EEG analysis by Scherg (for a review see Scherg 1990.),but the same approach canbe used in MEG studies as well (Scherg et al. 1989). The basic assumption of the model is that there are several dipolar sources which maintain their position and optionally also their orientation throughout the time interval of interest. However, the dipoles are allowed to change their amplitudes in order to produce a field distribution which matches the experimental values. We shall denote the measured and predicted data by the matrices ~ k and Bjlv respectively, where j = 1.....n indexes the measurement points and k = 1,...,m correspond to the time instants t k under consideration. We seek the minimum for the conventional least-squares error function S=
2
IIM-BII
(5)
F
w h e r e xl .....Xq are the u n k n o w n m o d e l parameters whereas I I. I I 2F denotes the square of the Frobenius norm:
t i ' a * t I2F = E
n
m
Ea~j
(6)
=
Tr(ara)
i=1 j=l We assume p dipoles located at r~, d = 1,...,p, and restrict our attention to the spherical conductor model so that any measurement is insensitive to radially oriented dipoles. There will be Pl fixed-orientation dipoles and p2 dipoles whose orientations vary as a function of time (p = Pl + P2). Therefore, we want to retrieve r = Pl + 2p2 dipole waveforms, one for each fixed-orientation dipole and two for each variable-orientation one. With these assumptions, the data predicted by the model can be written as: B (nx m) = G (nx 2p) ~1 ..... r~p R (2pxr) Q (rx m)
(7)
where the dimensions of each matrix appear as superscripts. Here G is a gain matrix composed of the unit dipole signals
Gj,2d_l= bj~d, e~) , Gj,2d bj (-~d, e-~), d = 1..... p, j= 1,...,u
(8)
where bj ( r ~ , e~ ) and bj ( r ~ , e~ ) are the signals that would be produced by unit dipoles at r~ pointing in the directions of the unit vectors e-~and e-7 of the spherical coordinate system, respectively. The matrix R contains the differentiation between fixed and variable-orientation dipoles:
MEG: A Tool for Functional Brain Imaging
"cos(~l)
0
......
99
...
sin(~l) 0 . . . . . . ... 0 cos(~2) 0 . . . . . . 0 sin(132) 0 . . . . . . R = ..,
...
0
...
0
......
0
...
...
...
0 0 0
0
cos(~pl) 0
0
sin(~pl ) 0
0
(9)
0
1(2p2 )
where I (2p2) is the 2p2 x 2p2 identity matrix. If all the orientations are varying, r -- 2p2 = 2p and R = I (2p). The fixed dipoles form angles ]3k, k = 1,..., Pl, with respect to e0. Finally, the first Pl rows of Q (see Eq. (7)) contain the amplitude time series of the fixed dipoles at tk, k = 1..... m, while the remaining 2p2 rows are the time series of the two components of each of the variable-orientation dipoles. Before minimizing Eq. (5), correct values must be chosen for Pl and P2. The proper choice can be made by comparing the singular values of the data matrix M with the estimated noise level of the measurements (de Munck 1990; Yin and Krishnaiah 1987; Wax and Kailath 1985). We shall, however, here omit the discussion of this important point and assume that Pl and P2 are known a priori. We shall focus on the case P2 = 0, whereby all dipoles have fixed orientations.
If all dipoles have fixed orientations, the angles 13k are usually included in the set of nonlinear parameters (Scherg 1990) and Q is found by the above linear approach. For this purpose, one just replaces G in Eqs. (10-12) G = G R. Recently it was shown (Mosher et al. 1990) that one could start with variable-orientation dipoles only and then find an optimal orientation for the fixed dipoles. In this procedure, Eq. (11) is first applied to find ~ for the r variable-orientation dipoles. Then, 2 x m submatrices ~k, k = 1,..., p, are formed, each containing two successive rows of ~ One then finds the best rank-one approximation forQk through the singular-value decomposition (Golub and van Loan 1989), viz., A
Qk
=
Ulk [ )~lk VIk T ]
(13) A
where ;~lk is the larger of the two singular values of Qk, while Ulk and Vlk are the corresponding left and right singular vectors, respectively. From Eq. (13) we can readily identify ul k = (cos 13k sin 13k)T. The term in brackets [. ] is the time series for the dipole.
Finding Dipole Positions The MUSIC a p p r o a c h
Linear Optimization If the model consists only of variable-orientation dipoles, R = ! and the data are given by B -- G (r'l .... r'p) Q. Assuming that the dipole locations r--d,d = 1..... p, are known, the remaining minimization problem Stain =
min II GQ - M II F2 Q
(10)
is linear with the solution a,
Q : G +M
(11)
where + denotes the pseudoinverse of a matrix (Golub and van Loan 1989). In fact, Smin =
[ I(I-GG+)MI
I F2 = I I P z M I
I F2
(12)
where P± is a projector to the left nullspace of G, see, for example (Golub and van Loan 1989). Minimizing over the nonlinear parameters r-~,d = 1..... p, Eq. (12) can be used to find Smin at each iteration step of the nonlinear optimization algorithm.
The remaining nonlinear optimization for r-~,d = 1,...,p, has been one of the central problems in multidipole modelling. This stems m a i n l y f r o m the fact that reasonable initial guesses are difficult to make. In addition, the calculation may get trapped into a local minim u m and Stain is never reached. So far the selection of the initial dipole positions has been based mainly on heuristic methods. Consequently, the results may strongly depend on the person who is in charge of the analysis. Especially if MEG is applied in clinical practice it is of vital importance to provide the means to remove these ambiguities. In the next two sections we shall concentrate on an elegant, recently proposed approach (Mosher et al. 1992) for searching the optimal locations of several simultaneously active current dipoles. This method is identical with the MUSIC algorithm (MUltiple Signal Classification), introduced earlier in another context (Schmidt 1986). The basic idea is to find the parameters giving rise to Stain in three steps: 1. Find the best possible projector, P± (Eq. (12)), regardless of G and R, assuming a certain combination of variable-orientation and fixed dipoles.
1O0
H~m61~inen
.R
3~
4 .,.~,.~
L
5 ,,..,,,,,./%~ 6 ,,-,,,,~P"%.~ 30 fT/cm I.___
100 msec 7--~/'~-
8.~,,~
9 , -J~J ~-~ 10,~'~" ~" z
40 mm
]
z
30 mm
25 t
Figure 3. Responses over the left hemisphere to electric stimulation of the contralateral thumb (interstimulus interval 2 s) m e a s u r e d with t h e Helsinki 2 4 - c h a n n e l gmdiometer (Ahonen et al, 1991) positioned over the left temporal and parietal lobe as shown in the inset. The lower and upper curves in each of the twelve locations show 8Bz/SX and 8Bz/Sy as a function of time. Here x, y, and z refer to the axes of a local coordinate system at each sensor location with z oriented perpendicular and x and y parallel to the curved bottom of the dewar. The passband is 0.05 - 90 Hz and 55 responses were averaged. The dashed lines show the data predicted by the two-dipole source model, with the dipoles located in the first and second somatosensory cortices as shown in figure 4. Modified from (Hari et al. 1992).
2. Find G and R for which I - (GR) (GR) + is as close as possible to the optimal projector. Equivalently, one can require that G R is orthonormal to Pz3. Find the dipole time series with the methods described in the previous section.
Fixed-orientation Dipole Scanning We assume next that all dipoles have fixed orientations: r = Pl- According to the recipe given in the previous section we first find the best orthogonal projector P± (minimizing Eq. 12) and then seek the set of dipole positions and orientations giving rise, as closely as possible, to this projector. We can accomplish this taskby evaluating a 'scanning function' at all feasible dipole positions. This function is inversely proportional to the probability that one of the dipoles is located at a given point on the basis of the information provided by our measurement. If our assumption about r fixed dipoles as a source is correct, we should find exactly r minima corresponding to the dipole locations we are searching for. The actual scanning function reads
ll~r~
j
~ 20
Qs,(t)
e~ ~s
100%
~
s o60
g(t) 0%
-50
-40
-30
-20
-10
0
Depth/rnm
Figure 4. (a) Contour plots of 1 / Sf (r) on two cross-sections of the head calculated from the data shown in figure 3. The approximate locations of the Rolandic and SyMan fissures are shown with gray lines. The two extrema correspond to the sources at SI (near the Rolandic fissure) and SII (near the Sylvian fissure). The origin of the z axis lies on the surface of the head; (b) 1 / Sf (r) plotted as a function of depth on radial lines passing through the two extrema seen in the z = -30 mm plot in (a); (c) The dipole waveforms and the goodness-of-fit calculated for the two optimal dipole positions deduced from (a) and (b). The stimulus onset time is indicated by a vertical arrow. Modified from (Hari et al. 1992).
Sf(
n'fin m TG 1 ( T P l G1 ( r-~m = m m TG1 ( T G1 ( r- m
(14)
where mk = (cos [~k sin ~k)T; G1 (r') is an n x 2 matrix containing the unit-dipole signals for a dipole at r-~,and P± is the optimal orthogonal projector given by P I = U IU IT
(15)
where U is an n x (m - r) matrix whose columns are the m - r left eigenvectors in the singular-value decomposition of the data matrix M, corresponding to its m - r smallest singular values. To find the dipole positions we can construct a contour plot of Sf ( r-~ or 1/Sf ( ~ and search for extrema. As an illustrative example we shall show here an application to the analysis of somatosensory evoked responses, see figures 3 and 4. If two fixed-orientation dipoles are assumed, the MUSIC scanning function peaks in this case at two locations which correspond
MEG: A Tool for Functional Brain Imaging
10t
However, the construction of the model requires identification of the shape of the cranium from MRI's which is not a straightforward task. With more complicated experimental paradigms it becomes possible to study interactions between multiple, spatially separated brain areas. For example, if it is expected that the visual cortex, the parietal cortex, and the auditory cortex could be activated in close sequence, measurements should be made over all these areas and processed together, in order to take into account the spatio-temporal overlaps correctly. To satisfy this goal it is necessary to use both a globally valid conductor model, instead of a locally fitted sphere, and a distributed source model, i.e., multiple dipoles or a continuous current distribution.
References
Figure 5. Superposition of the SI/SII dipole locations on 6-mm thick spin-echo (TR = 2500 ms, TE= 22 ms) MR images. In each case the image closest to the dipole location is selected. (a) and (b) The location of the SI source shown on transaxial and sagiYral images. (c) A coronal image showing both SI and SII sources. (d) Sagittal image showing the location of the SII source. Modified from Hari et al. 1992).
quite nicely to the locations of the first and second somatosensory cortices, see figure 5. In this case we were also able to confirm the MUSIC result simply by simple single-dipole fits to the data at time instants when one of the dipoles was silent (Hari et al. 1992).
Conclusions Because of the nonuniqueness of the neuromagnetic inverse problem, it is necessary to combine MEG data with supplementary information. Only the first steps in this direction have been taken. We need a system in which MRI gives structural information, PET and dynamic MRI provide metabolic and blood flow data, and EEG supplements MEG in obtaining data about signal processing in the brain. The improved accuracy of anatomic information and the combination of MEG with EEG will necessitate considerable improvements in modelling the conductivity of the head. The magnetic field can be accurately calculated from a realistically-shaped homogeneous head model; the process can be made computationally very effective.
Ahonen, A.I., H~im~il~iinen, M.S., Kajola, M.J., Knuutila, J.E.T., Lounasmaa, O.V., Simola, J.T., Tesche, C.D., and Vilkman, V.A. Multichannel SQUID systems for brain research. IEEE Trans. Magn., MAG, 1991, 27: 2786-2792. Arthur, R.M. and Geselowitz, D.B. Effect of inhomogeneities on the apparent location and magnitude of a cardiac current dipole source. IEEE Trans. Biomed. Eng., BME, 1970, 17: 141-146. Cuffin, B.N. and Cohen, D. Magnetic fields of a dipole in special volume conductor shapes. IEEE Trans. Biomed. Eng., BME, 1977, 24: 372-381. Geselowitz, D.B. On bioelectric potentials in an inhomogeneous volume conductor. Biophys. J., 1967~7:1--11. Golub, G.H. and van Loan, C.F. Matrix Computations. The Johns Hopkins University Press, Baltimore, second edition, 1989. Hari, R., Reinikainen, K., Kaukoranta, E., H~im~il~inen, M., I1moniemi, R., Penttinen, A., Salminen, J. and Teszner, D. Somatosensory evoked cerebral magnetic fields from SI and SII in man. Electroenceph. Clin. Neurophysiol., 1984, 57: 254-263. Hari, R. The neuromagnetic method in the study of the human auditory cortex. In F.Grandoli, M.Hoke, and G.L. Romani, editors, Auditory Evoked Magnetic Fields and Electric Potentials, Karger, Basel, 1990, 222-282. Hari, R. On brain's magnetic responses to sensory stimuli. J. Clin. Neurophysiol., 1991, 8: 157-169. Hari, R., Karhu, J., H~im~il~iinen, M., Knuutila, J., Sams, J., and Vilkman, V. Functional organization of the human first and second somatosensory cortex as revealed by neuromagnetic measurements. European J. Neurosci., 1992 (in press). Ham~ilainen, M.S., Hari, R., Ilmoniemi, R., Knuutila, J., and Lounasmaa, O.V. Magnetoencephalography - theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Mod. Phys., 1993, (in press). Ham~il~iinen M.S. and Ilmoniemi, R.J. Interpreting measured magnetic fields of the brain: Estimates of current distributions. Technical Report TKK-F-A559, Helsinki University of Technology, 1984. H~im~ilainen, M.S. and Sarvas, J. Realistic conductivity
102
geometry model of the human head for interpretation of neuromagnetic data. IEEE Trans. Biomed. Eng., 1989, 36: 165-171. Helmholtz, H. Ueber einige Gesetze der Vertheilung elektrischer StrOme in kOrperlichen Leitern, mit Anwendung auf die thierisch-elektrischen Versuche. Ann. Phys. Chem., 1853, 89: 211-233, 353-377. Hoenig, H.E., Daalmans, G.M., Bar, L, BOmmel, F., Paulus, A., Uhl, D., Weisse, H.J., Schneider, S., Seifert, H., Reichenberger, H., and Abraham-Fuchs, K. Multi channel dc SQUID sensor array for biomagnetic applications. IEEE Trans. Magn., MAG, 1991, 27: 2777-2785. Ilmoniemi, R.J., H~im~il~iinen, M.S. and Knuutila, J. The forward and inverse problems in the spherical model. In H. Weinberg, G. Stroink, and T. Katila, editors, Biomagnetism: Applications & Theory, Pergamon Press, New York, 1985, 278-282. Kajola, M., Ahlfors, S., Ehnholm, G.J., H~illstr6m, J., H~im~il~iinen, M.S., Ilmoniemi, R.J., Kiviranta, M., Knuutila, J., Lounasmaa, O.V., Tesche, C.D., and Vilkman, V. A 24channel magnetometer for brain research. In S.J. Williamson, M. Hoke, G. Stroink, and M. Kotani, editors, Advances in Biomagnetism, New York, Plenum, 1989, 673-676. Kaukoranta, E., Hari, R., H~im~il~iinen, M., and Huttunen, J. Cerebral magnetic fields evoked by peroneal nerve stimulation. Somatosensory Res., 1986, 3: 309-321. Kullmann, W.H., Jandt, K.D., Rehm, K., Schlitt, H.A., Dallas, W.J., and Smith, W.E. A linear estimation approach to biomagnetic imaging. In S.J. Williamson, M. Hoke, G. Stroink, and M. Kotani, editors, Advances in Biomagnetism, New York, Plenum, 1989, 571-574. Meijs, J.W.H., and Peters, M.J. The EEG and MEG, using a model of eccentric spheres to describe the head. IEEE Trans. Biomed. Eng., BME, 1987, 34: 913-920. Mosher, J.C., Lewis, P.S., and Leahy R. Multiple dipole modeling and localization from spatio-temporal MEG data. IEEE Trans. Biomed. Eng., 1992, BME-39: 541-557. Mosher, J.C., Lewis, P.S., Leahy, R. and Singh, M. Multiple dipole modeling of spatio-temporal MEG data. In A.F. Gmitro, P.S. Idell, and I.J. LaHaie, editors, Digital Image Synthesis and Inverse Optics, Proc. SPIE 1351,1990, 364-375. de Munck, J.C. The potential distribution in a layered spheroidal volume conductor. J. App1. Phys., 1988, 64: 464470. de Munck, J.C. A mathematical and physical interpretation of the electromagnetic field of the brain. PhD thesis, University of Amsterdam, 1989. de Munck, J.C. The estimation of time-varying dipoles on the basis of e v o k e d potentials. Electroenceph. Clin.
H~m~l~inen
Neurophysiol., 1990, 77: 156-160. de Munck, J.C. A linear discretization of the volume conductor boundary integral equation using analytically integrated elements. IEEE Tram: Biomed. Eng., in press 1992. Oostendorp, T.F. and van Oosterom, A. Source parameter estimation in inhomogeneous volume conductors of arbitrary shape. IEEE Trans. Biomed. Eng., BME, 1989, 36: 382-391. Ribary, U., Ioannides, A.A., Singh, K.D., Hasson, R., Bolton, J.P.R., Lado, F., Mogilner, A., and Llin'as, R. Magnetic field tomography of coherent thalamo-cortica140-Hz oscillations in humans. Proc. Natl. Acad. Sci. USA, 1991, 88: 1103711041. Romani, G.L., and Rossini, P. Neuromagnetic functional localization: Principles, state of the art, and perspectives. Brain Topography, 1988, 1: 5-19. Ryh~.nen, T., Sepp~i, H., Ilmoniemi, R., and Knuutila, J. SQUID magnetometers for low-frequency applications. J. Low Temp. Phys., 1989, 76: 287-386. Sarvas, J. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. Phys. Med. Biol., 1987, 32: 11-22. Scherg, M., Hari, R., and H~im~il~iinen, M. Frequency-specific sources of the auditory N19-P30-P50 response detected by a multiple source analysis of evoked magnetic fields and potentials. In S.J. Williamson, M. Hoke, G. Stroink, and M. Kotani, editors, Advances in Biomagnetism, New York, Plenum, 1989, 97-100. Scherg, M. Fundamentals of dipole source potential analysis. Karger, Basel, 1990, 40-69. Schmidt, R.O. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propagat., AP, 1986, 34: 276-280. Urankar, L. Common compact analytical formulas for computation of geometry integrals on a basic Cartesian subdomain in boundary and volume integral methods. Eng. Anal. Boundary Elements, 1990, 7: 124-129. Wax, M., and Kailath, T. Detection of signals by information theoretic criteria. IEEE Trans. Acoust., Speech and Signal Processing, ASSP, 1985, 33: 387-392. Williamson, S.J., and Kaufman, L. Biomagnetism. J. Magn. Magn. Mat., 1981, 22: 129-201. Williamson, S.J., Lii, Z.-L., Karron, D. and Kaufman, L Advantages and limitations of magnetic source imaging. Brain Topography, 1991, 4: 169-180. Yin, Y.Q. and Krishnaiah, P.R. On some nonparametric methods for detection of the number of signals. IEEE Trans. Acoust., Speech and Signal Processing, ASSP, 1987, 35:15331538.
95
Magnetoencephalography: A Tool For Functional Brain Imaging Matti S. H~m~Jl~inen* Summary: At present, one of the most promising windowsto the functionalorganizationof the human brain is magnetoencephalography(MEG). Bymappingthe magneticfielddistributionoutsidethe head the sites of neuraleventscan be locatedwith an accuracyof a few millimetersand the temporalevolutionofthe activationcanbe tracedwitha millisecondresolution.Thispaper reviewssomeforwardfieldcalculationapproachessuitable forthe interpretationofthe brain's electromagneticsignals.Inversemodellingwithmultipledipolesis describedin detail.An exampleof the analysis ofthe sornatosensoryevoked-responsesillustratesthepotentialofmultiplesignalclassification(MUSIC)algorithmin findingoptimaldipolepositions. Keywords: Magnetoencephalography;Functionalbrain imaging;Forwardproblem; Currentsource estimation.
Introduction Modern imaging techniques provide high-resolution windows through which the structure and metabolism of the h u m a n brain can be studied in vivo. Computer-assisted X-ray tomography (CT) and magnetic-resonance imaging (MRI) produce high-quality but static maps of the brain's anatomy. Functional information about the brain on a time scale of minutes can be obtained with regional blood flow (RCBF) and positron-emission tomography (PET) measurements, based on changes in cell metabolism as indicated by radioactive tracers. There are only two totally noninvasive windows to study the actual brain function on a millisecond scale: electroencephalography (EEG) and magnetoencephalography (MEG). EEG, employed for about 60 years, involves the measurement of the electric potential distribution by means of electrodes attached on the scalp. In MEG, the minute (50 - 500fT) magnetic fields, produced by the same brain events as in EEG, are detected outside the head with SQUID (Superconducting QUantum Interference Device) magnetometers. The processes which can be studied by means of MEG *LowTemperatureLaboratory,HelsinkiUniversityof Technology, Espoo,Finland. Acceptedfor publication: June23, 1992. I thank Prof. Riitta Hari, Prof. Olli V. Lounasmaa,and Dr. Juha Simola for valuable comments on this manuscript. This work was supported by Finnish National Fund for Researchand Development (SITRA)and InstrumentariumCorp. Correspondenceand reprint requests shouldbe addressedto Matti S. Hamalainen,Low Temperature Laboratory,HelsinkiUniversityof Technology,02150Espoo,Finland. Copyright© 1992Human SciencesPress, Inc.
may be divided into two main categories: evoked responses and spontaneous brain activity. In an evoked-field experiment, sensory stimuli are presented to the subject and the ensuing magnetic field is recorded outside the head at several locations; at each site the stimulus is repeated 50 - 200 times in order to remove incoherent noise by signal averaging. Spontaneous phenomena, both normal and abnormal, resulting from localized events in the brain, e.g., the 10 Hz alpha-rhythm and epileptic activity, can also be investigated provided that the signal-to-noise ratio is sufficient, The first MEG measurements were made with induction coil magnetometers in the late 1960's. At present, multichannel SQUID devices of considerable sophistication are used by several research groups around the world. These instruments include the Helsinki 24-channel system (Kajola et al. 1989; Ahonen et al. 1991), the Siemens 37-channel 'KrenikonTM' (Hoenig et al. 1991), and the 37-channel 'Magnes TM' of Biomagnetic Technologies, Inc. Although the present instruments allow several types of experiments to be carried out, it has become evident that, in order to study complicated cognitive processes and s p o n t a n e o u s activity in h e a l t h y subjects and patients, it is necessary to record the magnetic field over the whole scalp simultaneously. The first pilot evokedresponse measurements with a 122-channel 'Neuromag122' system have already been carried out in our laboratory. Rather than attempting to cover all aspects of MEG recordings and their interpretation, we will focus on the problem of neuromagnetic source modelling, emphasizing the new requirements due to simultaneous full-head
96
H6m~l~inen
proble~ Forward
Maxwetl*s equations (~
Inverse problem
Figure 1. Schematic illustration of the forward and inverse problems.
coverage. We shall elaborate on the different available models for forward field computation and for the source itself. For in-depth treatises on neuromagnetic measurements and their interpretation an interested reader may wish to consult some of the available review articles (Williamson and Kaufman 1981; Sarvas 1987; Romani and Rossini 1988; Ryh~inen et al. 1989; Hari 1991; Williamson et al. 1991; H~im~il~iinenet al. 1993).
Source Modelling The neuromagnetic inverse problem is to estimate the cerebral current sources underlying a measured distribution of the magnetic field (see figure 1). We shall limit our discussion to this task and will not elaborate on the physiological interpretation, although we shall keep in mind that neuromagnetic fields arise in the brain, which imposes constraints on the sources that the measurements might reveal. Contrary to CT and MRI, the inverse problem of MEG does not have a unique solution (Helrnholtz 1853), even if the magnetic field were precisely known everywhere outside the head. Inclusion of EEG information does not help: we must still restrict the viable source current distributions by choosing a suitable source model. In fact, one can find current distributions that are either magnetically silent (MEG vanishes), electrically silent (EEG vanishes), or both. A simple example of a magnetically silent source that produces an electric field is a radial dipole in a spherically symmetric conductor. The opposite is a current loop, which is electrically silent but produces a magnetic field. Because of the nonuniqueness of the inverse problem it is necessary to confine to a restricted class of source currents. We shall here discuss the problems associated with modelling the neural activity as a collection of point sources, current dipoles. An alternative approach of using distributed source m o d e l l i n g was also i n t r o d u c e d at an early stage (Ham~il~iinen and Ilmoniemi 1984; Ilmoniemi et al. 1985) but, due to technical limitations, it has become useful only after the introduction of the multichannel magnetometers (Kullmann et al. 1989; Ribary et al. 1991).
Forward solutions The solution of the inverse problem necessarily involves the forward field and potential calculations. We shall, therefore, next give a synopsis of the results involved in the computation of the electromagnetic field ensuing from neural current sources. We shall also point out some characteristic differences between MEG and EEG, which emerge from the theoretical results. From a known distribution of currents it is, in principle, straightforward to compute the corresponding distribution of the magnetic field and electric potential using Maxwell's equations. To accomplish the calculation it is necessary to know the spatial distribution of the electrical conductivity, the dielectric constant, and the magnetic permeability. However, for MEG and EEG the calculation is simplified because we can omit the time dependent terms from Maxwell's equations: the quasistatic approximation is valid. Furthermore, we can assume that the magnetic permeabilities of the media in question are equal to that of free space. The remaining problem of field computation is illustrated in figure 1. We start from a known neural current distribution, the primary current J--~,which gives rise to a passive ohmic current, commonly called the volume current, P , extending throughout the head. The volume current is directly related to the electric t~otential gradient by the Ohm s law: jv ( r-~ = - o( ~ vV( r--J,where is the conductivity and V the electric potential; EEG involves the measurement of V on the surface of the head. It is w o r t h n o t i n g that the c o n d u c t i v i t y m a y be anisotropic. Computation of the electric potential is then numerically more complex than in the isotropic case as will be discussed later in this article. Once we know both J--~,and J-~, the magnetic field can be calculated from the Biot-Savart's law:
(1)
MEG: A Tool for Functional Brain Imaging
97
Measurernsnt I
Figure 2. The choice of the origin in the sphere model should be based on the local curvature of the skull's inner surface. The two measurement areas shown require different origins, O1 and 02. The corresponding radii of curvature of the two areas are denoted by R1 and R2, respectively. Typical sources that we would expect to locate with _~easur~ments over the areas shown are denoted by Q1 and Q2. (right) A piecewise homogneous conductor. Each region Gi has uniform conductivity ~i. Unit vectors normal to the surfaces are denoted by nij.
The Sphere Model If the conductor has__a simple symmetrical shape, the computation the V and B can become remarkably simple: we may obtain quickly converging series representations or even an analytic expression. The sphere model is often used in MEG and EEG forward field calculations: the head is assumed to consist of several spherically symmetric layers. It turns out that the magnetic field remains the same, irrespective of the radii and conductivities of the shells. As discussed by de Munck (1989), this is the case even in an anisotropic conductor if there is a difference between the radial and the two tangential conducfivities only. One can derive an analytical formulae for B allowing fast evaluation. The sphere model can be readily applied once we have selected a suitable center of symmetry. Since the skull is a poor conductor, major contribution to the magnetic field comes from the currents flowing within the cranial space. The sphere model should, therefore, be fitted to the inner curvature of the skull, as shown in figure 2. This task is straightforward if MR/images from the subject are available. Otherwise, one must rely on an approximation from the outer shape of the head. From the sphere model one can predict that radial sources are not seen in MEG and that deep sources are strongly attenuated. This means that MEG selectively detects the activity of fissural cortex. Fortunately, many interesting areas, such as the primary auditory, motor, somatosensory, and visual cortices, are located in fissures. Sources in them are~ therefore, easily detectable by MEG. Furthermore, the activity of the convexial cortex,
visible in EEG, is filtered out, whereby the signal-to-noise ratio is improved. The electric potential on the surface of a spherically symmetric conductor is affected both by the conductivity profile and anisotropies in conductivity. For example, if a poorly conducting layer (skull) is placed between two relatively good conductors (brain and scalp), the potential pattern caused by a current source inside will be more widespread than in a homogeneous sphere. Another difference is that both radial and deep sources contribute significantly to V. The computation of V requires the evaluation of a series expansion in spherical harmonics (Arthur and Geselowitz 1970; de Munck 1988), which is more time consuming than the computation of B from simple analytical expressions.
Modelling Conductors of Arbitrary Shape Other simple conductor shapes that can be handled analytically are spheroids (de Munck 1988; Cuffin and Cohen 1977) and eccentric spheres (Meijs and Peters 1987). However, if the shape of the conductor is not restricted we must resort to more general numerical methods for computing V and B. In the quasistatic approximation, the divergence of the total current density ~r--) = J-~ (r--) + j g ( r ' ) vanishes: V. J-~ r--) : O, which leads directly to the Poisson's equation for the electric potential: v. (
vv) = v. pp
(2)
With cy and J~ known, we can solve this differential equation numerically by means of the finite element method (FEM). This approach is very general: we can have any conductivity distribution with anisotropies. We shall focus on another possibility, the boundaryelement method. The head is then assumed to consist of several regions with uniform and isotropic conductivities, see figure 2. Using Green's identity it is now possible to convert the differential equation for V into an integral equation involving V on the boundaries of regions of different conductivity only. For a single homogeneous, arbitrarily shaped conductor with conductivity cy,surrounded by an insulator, the equation reads (Geselowitz 1967):
V ( ~ = 2Vo(r-} + G
(3)
where r' is a point on the bounding surface S, V0( ~ is the potential if the conductor had an infinite extent, and dr2 r-( r-+) is the solid angle subtended at r---boya surface element at r-9'.
98
H6m~l~inen
For a solution of this equation, S is approximated by a net of triangles (Oostendorp and van Oosterom 1989; H~im~ilainen and Sarvas 1989). We can then either assume that V is constant within one triangle or use higher-order approximations to yield more accurate results with smaller computational burden (Urankar 1990; de Munck 1992). The discretization results in a set of linear algebraic equations which can be solved for the unknown potential. Once V is known on S, we can use the Geselowitz' formula, derived from Eq. (1):
N0
v ( ~ ' ) v"
1
tr~_r-+,t
x d ~ r--~')
(4) to compute the magnetic field. In one series of simulations (H~im~il~iinenand Sarvas 1989) comparing the sphere model with more accurate, realistically shaped models, it was found that, except for frontotemporal sources, the sphere provides a sufficient accuracy for practical purposes. When the sphere model is imprecise, a brain-shaped homogeneous model, which assumes that the skull is a perfect insulator, is an excellent approximation. It must be noted that if one likes to predict the EEG potential distribution with the same accuracy as that provided by the homogeneous model in MEG, one must apply a multilayer realistically shaped head model with known conductivities and shapes for the brain, the skull, and the scalp. Much less information is thus needed for the analysis of MEG data.
Multidipole Models Formulation of the Problem
An obvious generalization of the popular singledipole model is an assumption of multiple spatially separated dipolar sources. If the distance between the individual dipoles is sufficiently large ( >3 cm) and their orientations are favorable, the field patterns may show only minor overlap and they can be fitted individually using the single-dipole model. An example of this approach is the separation of activities from the first and s e c o n d s o m a t o s e n s o r y cortices (Hari et al. 1984; Kaukoranta et al. 1986; Hari 1990). Similarly, if the temporal behavior of the dipoles differ, it is often possible to recognize each source separately. However, when the sources overlap both temporally and spatially, we must resort to the multidipole approach to obtain correct results. An effective approach to this modelling problem is to take into account the spatiotemporal course of the signals as a whole, instead of considering each time sample separately. This method was first
applied to EEG analysis by Scherg (for a review see Scherg 1990.),but the same approach canbe used in MEG studies as well (Scherg et al. 1989). The basic assumption of the model is that there are several dipolar sources which maintain their position and optionally also their orientation throughout the time interval of interest. However, the dipoles are allowed to change their amplitudes in order to produce a field distribution which matches the experimental values. We shall denote the measured and predicted data by the matrices ~ k and Bjlv respectively, where j = 1.....n indexes the measurement points and k = 1,...,m correspond to the time instants t k under consideration. We seek the minimum for the conventional least-squares error function S=
2
IIM-BII
(5)
F
w h e r e xl .....Xq are the u n k n o w n m o d e l parameters whereas I I. I I 2F denotes the square of the Frobenius norm:
t i ' a * t I2F = E
n
m
Ea~j
(6)
=
Tr(ara)
i=1 j=l We assume p dipoles located at r~, d = 1,...,p, and restrict our attention to the spherical conductor model so that any measurement is insensitive to radially oriented dipoles. There will be Pl fixed-orientation dipoles and p2 dipoles whose orientations vary as a function of time (p = Pl + P2). Therefore, we want to retrieve r = Pl + 2p2 dipole waveforms, one for each fixed-orientation dipole and two for each variable-orientation one. With these assumptions, the data predicted by the model can be written as: B (nx m) = G (nx 2p) ~1 ..... r~p R (2pxr) Q (rx m)
(7)
where the dimensions of each matrix appear as superscripts. Here G is a gain matrix composed of the unit dipole signals
Gj,2d_l= bj~d, e~) , Gj,2d bj (-~d, e-~), d = 1..... p, j= 1,...,u
(8)
where bj ( r ~ , e~ ) and bj ( r ~ , e~ ) are the signals that would be produced by unit dipoles at r~ pointing in the directions of the unit vectors e-~and e-7 of the spherical coordinate system, respectively. The matrix R contains the differentiation between fixed and variable-orientation dipoles:
MEG: A Tool for Functional Brain Imaging
"cos(~l)
0
......
99
...
sin(~l) 0 . . . . . . ... 0 cos(~2) 0 . . . . . . 0 sin(132) 0 . . . . . . R = ..,
...
0
...
0
......
0
...
...
...
0 0 0
0
cos(~pl) 0
0
sin(~pl ) 0
0
(9)
0
1(2p2 )
where I (2p2) is the 2p2 x 2p2 identity matrix. If all the orientations are varying, r -- 2p2 = 2p and R = I (2p). The fixed dipoles form angles ]3k, k = 1,..., Pl, with respect to e0. Finally, the first Pl rows of Q (see Eq. (7)) contain the amplitude time series of the fixed dipoles at tk, k = 1..... m, while the remaining 2p2 rows are the time series of the two components of each of the variable-orientation dipoles. Before minimizing Eq. (5), correct values must be chosen for Pl and P2. The proper choice can be made by comparing the singular values of the data matrix M with the estimated noise level of the measurements (de Munck 1990; Yin and Krishnaiah 1987; Wax and Kailath 1985). We shall, however, here omit the discussion of this important point and assume that Pl and P2 are known a priori. We shall focus on the case P2 = 0, whereby all dipoles have fixed orientations.
If all dipoles have fixed orientations, the angles 13k are usually included in the set of nonlinear parameters (Scherg 1990) and Q is found by the above linear approach. For this purpose, one just replaces G in Eqs. (10-12) G = G R. Recently it was shown (Mosher et al. 1990) that one could start with variable-orientation dipoles only and then find an optimal orientation for the fixed dipoles. In this procedure, Eq. (11) is first applied to find ~ for the r variable-orientation dipoles. Then, 2 x m submatrices ~k, k = 1,..., p, are formed, each containing two successive rows of ~ One then finds the best rank-one approximation forQk through the singular-value decomposition (Golub and van Loan 1989), viz., A
Qk
=
Ulk [ )~lk VIk T ]
(13) A
where ;~lk is the larger of the two singular values of Qk, while Ulk and Vlk are the corresponding left and right singular vectors, respectively. From Eq. (13) we can readily identify ul k = (cos 13k sin 13k)T. The term in brackets [. ] is the time series for the dipole.
Finding Dipole Positions The MUSIC a p p r o a c h
Linear Optimization If the model consists only of variable-orientation dipoles, R = ! and the data are given by B -- G (r'l .... r'p) Q. Assuming that the dipole locations r--d,d = 1..... p, are known, the remaining minimization problem Stain =
min II GQ - M II F2 Q
(10)
is linear with the solution a,
Q : G +M
(11)
where + denotes the pseudoinverse of a matrix (Golub and van Loan 1989). In fact, Smin =
[ I(I-GG+)MI
I F2 = I I P z M I
I F2
(12)
where P± is a projector to the left nullspace of G, see, for example (Golub and van Loan 1989). Minimizing over the nonlinear parameters r-~,d = 1..... p, Eq. (12) can be used to find Smin at each iteration step of the nonlinear optimization algorithm.
The remaining nonlinear optimization for r-~,d = 1,...,p, has been one of the central problems in multidipole modelling. This stems m a i n l y f r o m the fact that reasonable initial guesses are difficult to make. In addition, the calculation may get trapped into a local minim u m and Stain is never reached. So far the selection of the initial dipole positions has been based mainly on heuristic methods. Consequently, the results may strongly depend on the person who is in charge of the analysis. Especially if MEG is applied in clinical practice it is of vital importance to provide the means to remove these ambiguities. In the next two sections we shall concentrate on an elegant, recently proposed approach (Mosher et al. 1992) for searching the optimal locations of several simultaneously active current dipoles. This method is identical with the MUSIC algorithm (MUltiple Signal Classification), introduced earlier in another context (Schmidt 1986). The basic idea is to find the parameters giving rise to Stain in three steps: 1. Find the best possible projector, P± (Eq. (12)), regardless of G and R, assuming a certain combination of variable-orientation and fixed dipoles.
1O0
H~m61~inen
.R
3~
4 .,.~,.~
L
5 ,,..,,,,,./%~ 6 ,,-,,,,~P"%.~ 30 fT/cm I.___
100 msec 7--~/'~-
8.~,,~
9 , -J~J ~-~ 10,~'~" ~" z
40 mm
]
z
30 mm
25 t
Figure 3. Responses over the left hemisphere to electric stimulation of the contralateral thumb (interstimulus interval 2 s) m e a s u r e d with t h e Helsinki 2 4 - c h a n n e l gmdiometer (Ahonen et al, 1991) positioned over the left temporal and parietal lobe as shown in the inset. The lower and upper curves in each of the twelve locations show 8Bz/SX and 8Bz/Sy as a function of time. Here x, y, and z refer to the axes of a local coordinate system at each sensor location with z oriented perpendicular and x and y parallel to the curved bottom of the dewar. The passband is 0.05 - 90 Hz and 55 responses were averaged. The dashed lines show the data predicted by the two-dipole source model, with the dipoles located in the first and second somatosensory cortices as shown in figure 4. Modified from (Hari et al. 1992).
2. Find G and R for which I - (GR) (GR) + is as close as possible to the optimal projector. Equivalently, one can require that G R is orthonormal to Pz3. Find the dipole time series with the methods described in the previous section.
Fixed-orientation Dipole Scanning We assume next that all dipoles have fixed orientations: r = Pl- According to the recipe given in the previous section we first find the best orthogonal projector P± (minimizing Eq. 12) and then seek the set of dipole positions and orientations giving rise, as closely as possible, to this projector. We can accomplish this taskby evaluating a 'scanning function' at all feasible dipole positions. This function is inversely proportional to the probability that one of the dipoles is located at a given point on the basis of the information provided by our measurement. If our assumption about r fixed dipoles as a source is correct, we should find exactly r minima corresponding to the dipole locations we are searching for. The actual scanning function reads
ll~r~
j
~ 20
Qs,(t)
e~ ~s
100%
~
s o60
g(t) 0%
-50
-40
-30
-20
-10
0
Depth/rnm
Figure 4. (a) Contour plots of 1 / Sf (r) on two cross-sections of the head calculated from the data shown in figure 3. The approximate locations of the Rolandic and SyMan fissures are shown with gray lines. The two extrema correspond to the sources at SI (near the Rolandic fissure) and SII (near the Sylvian fissure). The origin of the z axis lies on the surface of the head; (b) 1 / Sf (r) plotted as a function of depth on radial lines passing through the two extrema seen in the z = -30 mm plot in (a); (c) The dipole waveforms and the goodness-of-fit calculated for the two optimal dipole positions deduced from (a) and (b). The stimulus onset time is indicated by a vertical arrow. Modified from (Hari et al. 1992).
Sf(
n'fin m TG 1 ( T P l G1 ( r-~m = m m TG1 ( T G1 ( r- m
(14)
where mk = (cos [~k sin ~k)T; G1 (r') is an n x 2 matrix containing the unit-dipole signals for a dipole at r-~,and P± is the optimal orthogonal projector given by P I = U IU IT
(15)
where U is an n x (m - r) matrix whose columns are the m - r left eigenvectors in the singular-value decomposition of the data matrix M, corresponding to its m - r smallest singular values. To find the dipole positions we can construct a contour plot of Sf ( r-~ or 1/Sf ( ~ and search for extrema. As an illustrative example we shall show here an application to the analysis of somatosensory evoked responses, see figures 3 and 4. If two fixed-orientation dipoles are assumed, the MUSIC scanning function peaks in this case at two locations which correspond
MEG: A Tool for Functional Brain Imaging
10t
However, the construction of the model requires identification of the shape of the cranium from MRI's which is not a straightforward task. With more complicated experimental paradigms it becomes possible to study interactions between multiple, spatially separated brain areas. For example, if it is expected that the visual cortex, the parietal cortex, and the auditory cortex could be activated in close sequence, measurements should be made over all these areas and processed together, in order to take into account the spatio-temporal overlaps correctly. To satisfy this goal it is necessary to use both a globally valid conductor model, instead of a locally fitted sphere, and a distributed source model, i.e., multiple dipoles or a continuous current distribution.
References
Figure 5. Superposition of the SI/SII dipole locations on 6-mm thick spin-echo (TR = 2500 ms, TE= 22 ms) MR images. In each case the image closest to the dipole location is selected. (a) and (b) The location of the SI source shown on transaxial and sagiYral images. (c) A coronal image showing both SI and SII sources. (d) Sagittal image showing the location of the SII source. Modified from Hari et al. 1992).
quite nicely to the locations of the first and second somatosensory cortices, see figure 5. In this case we were also able to confirm the MUSIC result simply by simple single-dipole fits to the data at time instants when one of the dipoles was silent (Hari et al. 1992).
Conclusions Because of the nonuniqueness of the neuromagnetic inverse problem, it is necessary to combine MEG data with supplementary information. Only the first steps in this direction have been taken. We need a system in which MRI gives structural information, PET and dynamic MRI provide metabolic and blood flow data, and EEG supplements MEG in obtaining data about signal processing in the brain. The improved accuracy of anatomic information and the combination of MEG with EEG will necessitate considerable improvements in modelling the conductivity of the head. The magnetic field can be accurately calculated from a realistically-shaped homogeneous head model; the process can be made computationally very effective.
Ahonen, A.I., H~im~il~iinen, M.S., Kajola, M.J., Knuutila, J.E.T., Lounasmaa, O.V., Simola, J.T., Tesche, C.D., and Vilkman, V.A. Multichannel SQUID systems for brain research. IEEE Trans. Magn., MAG, 1991, 27: 2786-2792. Arthur, R.M. and Geselowitz, D.B. Effect of inhomogeneities on the apparent location and magnitude of a cardiac current dipole source. IEEE Trans. Biomed. Eng., BME, 1970, 17: 141-146. Cuffin, B.N. and Cohen, D. Magnetic fields of a dipole in special volume conductor shapes. IEEE Trans. Biomed. Eng., BME, 1977, 24: 372-381. Geselowitz, D.B. On bioelectric potentials in an inhomogeneous volume conductor. Biophys. J., 1967~7:1--11. Golub, G.H. and van Loan, C.F. Matrix Computations. The Johns Hopkins University Press, Baltimore, second edition, 1989. Hari, R., Reinikainen, K., Kaukoranta, E., H~im~il~inen, M., I1moniemi, R., Penttinen, A., Salminen, J. and Teszner, D. Somatosensory evoked cerebral magnetic fields from SI and SII in man. Electroenceph. Clin. Neurophysiol., 1984, 57: 254-263. Hari, R. The neuromagnetic method in the study of the human auditory cortex. In F.Grandoli, M.Hoke, and G.L. Romani, editors, Auditory Evoked Magnetic Fields and Electric Potentials, Karger, Basel, 1990, 222-282. Hari, R. On brain's magnetic responses to sensory stimuli. J. Clin. Neurophysiol., 1991, 8: 157-169. Hari, R., Karhu, J., H~im~il~iinen, M., Knuutila, J., Sams, J., and Vilkman, V. Functional organization of the human first and second somatosensory cortex as revealed by neuromagnetic measurements. European J. Neurosci., 1992 (in press). Ham~ilainen, M.S., Hari, R., Ilmoniemi, R., Knuutila, J., and Lounasmaa, O.V. Magnetoencephalography - theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Mod. Phys., 1993, (in press). Ham~il~iinen M.S. and Ilmoniemi, R.J. Interpreting measured magnetic fields of the brain: Estimates of current distributions. Technical Report TKK-F-A559, Helsinki University of Technology, 1984. H~im~ilainen, M.S. and Sarvas, J. Realistic conductivity
102
geometry model of the human head for interpretation of neuromagnetic data. IEEE Trans. Biomed. Eng., 1989, 36: 165-171. Helmholtz, H. Ueber einige Gesetze der Vertheilung elektrischer StrOme in kOrperlichen Leitern, mit Anwendung auf die thierisch-elektrischen Versuche. Ann. Phys. Chem., 1853, 89: 211-233, 353-377. Hoenig, H.E., Daalmans, G.M., Bar, L, BOmmel, F., Paulus, A., Uhl, D., Weisse, H.J., Schneider, S., Seifert, H., Reichenberger, H., and Abraham-Fuchs, K. Multi channel dc SQUID sensor array for biomagnetic applications. IEEE Trans. Magn., MAG, 1991, 27: 2777-2785. Ilmoniemi, R.J., H~im~il~iinen, M.S. and Knuutila, J. The forward and inverse problems in the spherical model. In H. Weinberg, G. Stroink, and T. Katila, editors, Biomagnetism: Applications & Theory, Pergamon Press, New York, 1985, 278-282. Kajola, M., Ahlfors, S., Ehnholm, G.J., H~illstr6m, J., H~im~il~iinen, M.S., Ilmoniemi, R.J., Kiviranta, M., Knuutila, J., Lounasmaa, O.V., Tesche, C.D., and Vilkman, V. A 24channel magnetometer for brain research. In S.J. Williamson, M. Hoke, G. Stroink, and M. Kotani, editors, Advances in Biomagnetism, New York, Plenum, 1989, 673-676. Kaukoranta, E., Hari, R., H~im~il~iinen, M., and Huttunen, J. Cerebral magnetic fields evoked by peroneal nerve stimulation. Somatosensory Res., 1986, 3: 309-321. Kullmann, W.H., Jandt, K.D., Rehm, K., Schlitt, H.A., Dallas, W.J., and Smith, W.E. A linear estimation approach to biomagnetic imaging. In S.J. Williamson, M. Hoke, G. Stroink, and M. Kotani, editors, Advances in Biomagnetism, New York, Plenum, 1989, 571-574. Meijs, J.W.H., and Peters, M.J. The EEG and MEG, using a model of eccentric spheres to describe the head. IEEE Trans. Biomed. Eng., BME, 1987, 34: 913-920. Mosher, J.C., Lewis, P.S., and Leahy R. Multiple dipole modeling and localization from spatio-temporal MEG data. IEEE Trans. Biomed. Eng., 1992, BME-39: 541-557. Mosher, J.C., Lewis, P.S., Leahy, R. and Singh, M. Multiple dipole modeling of spatio-temporal MEG data. In A.F. Gmitro, P.S. Idell, and I.J. LaHaie, editors, Digital Image Synthesis and Inverse Optics, Proc. SPIE 1351,1990, 364-375. de Munck, J.C. The potential distribution in a layered spheroidal volume conductor. J. App1. Phys., 1988, 64: 464470. de Munck, J.C. A mathematical and physical interpretation of the electromagnetic field of the brain. PhD thesis, University of Amsterdam, 1989. de Munck, J.C. The estimation of time-varying dipoles on the basis of e v o k e d potentials. Electroenceph. Clin.
H~m~l~inen
Neurophysiol., 1990, 77: 156-160. de Munck, J.C. A linear discretization of the volume conductor boundary integral equation using analytically integrated elements. IEEE Tram: Biomed. Eng., in press 1992. Oostendorp, T.F. and van Oosterom, A. Source parameter estimation in inhomogeneous volume conductors of arbitrary shape. IEEE Trans. Biomed. Eng., BME, 1989, 36: 382-391. Ribary, U., Ioannides, A.A., Singh, K.D., Hasson, R., Bolton, J.P.R., Lado, F., Mogilner, A., and Llin'as, R. Magnetic field tomography of coherent thalamo-cortica140-Hz oscillations in humans. Proc. Natl. Acad. Sci. USA, 1991, 88: 1103711041. Romani, G.L., and Rossini, P. Neuromagnetic functional localization: Principles, state of the art, and perspectives. Brain Topography, 1988, 1: 5-19. Ryh~.nen, T., Sepp~i, H., Ilmoniemi, R., and Knuutila, J. SQUID magnetometers for low-frequency applications. J. Low Temp. Phys., 1989, 76: 287-386. Sarvas, J. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. Phys. Med. Biol., 1987, 32: 11-22. Scherg, M., Hari, R., and H~im~il~iinen, M. Frequency-specific sources of the auditory N19-P30-P50 response detected by a multiple source analysis of evoked magnetic fields and potentials. In S.J. Williamson, M. Hoke, G. Stroink, and M. Kotani, editors, Advances in Biomagnetism, New York, Plenum, 1989, 97-100. Scherg, M. Fundamentals of dipole source potential analysis. Karger, Basel, 1990, 40-69. Schmidt, R.O. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propagat., AP, 1986, 34: 276-280. Urankar, L. Common compact analytical formulas for computation of geometry integrals on a basic Cartesian subdomain in boundary and volume integral methods. Eng. Anal. Boundary Elements, 1990, 7: 124-129. Wax, M., and Kailath, T. Detection of signals by information theoretic criteria. IEEE Trans. Acoust., Speech and Signal Processing, ASSP, 1985, 33: 387-392. Williamson, S.J., and Kaufman, L. Biomagnetism. J. Magn. Magn. Mat., 1981, 22: 129-201. Williamson, S.J., Lii, Z.-L., Karron, D. and Kaufman, L Advantages and limitations of magnetic source imaging. Brain Topography, 1991, 4: 169-180. Yin, Y.Q. and Krishnaiah, P.R. On some nonparametric methods for detection of the number of signals. IEEE Trans. Acoust., Speech and Signal Processing, ASSP, 1987, 35:15331538.
