*Nicolet Biomedical Inc., Madison, WI, U.S.A.. ^Department of ... of Wisconsin - Madison, Madison, WI, U.S.A.. #Department ... (1992) and Spencer et al. (1992).
Brain Topography, Volume 9, Number 1, t996
39
A Spatial Filtering Technique to Detect and Localize Multiple Sources in the Brain W. van Drongelen*, M. Yuchtman*, B. D. Van Veen ^, and A.C. van Huffelen #
Summary: An algorithm for localization of electromagnetic activity in the central nervous system is explored. This algorithm generates a neural activity index map within the brain by passing surface recordings through a set of spatial filters. The covariance matrix of the surface recordings is used to optimize the spatial filters' responses. This approach is studied in simulated situations and in real data. The simulations show the method's capability to detect areas of activity without prior knowledge of the number of sources. The resolving power of the method increases with number of electrodes and signal-to-noise ratio, and it decreases with depth. The analysis of the electrophysiologicaldata indicates that the method can distinguish simultaneously active areas in a realistic fashion. The analyzed recordings are bilateral median SEP responses, an epoch of spike activity showing several active regions and a recording with eye movement superimposed on spike activity. The method and the results are discussed in relation to current localization techniques. Key words: Spatial filtering; Source localization;Multiple sources; Spike localization;N20.
Introduction Localization of neural activity in the brain can generate insight into basic and clinical aspects of neural processing. Finding the anatomical correlate of epileptiform activity in the EEG is important clinically (Wong and Weinberg 1988; Ebersole 1991; van der Meij 1992; Scherg and Ebersole 1994). Especially for surgical candidates, the localization of epileptic discharges requires b o t h high reliability and resolution. For the reasons described, m a n y theoretical studies examine the relationship between source location and recorded electromagnetic activity (e.g., Scherg and v o n C r a m o n 1985; Fender 1987; K a v a n a g h et al. 1987; Meijs 1988; N u n e z 1988; H a m a l a i n e n a n d Sarvas 1989;de Munck 1989; Salu et al. 1990; Scherg 1990, 1992; Yan et al. 1991; Mosher et al. 1992; v a n Veen et al. 1992; Dale and Sereno 1993; Ioannides et al. 1993; Roth et al. 1993; Amir 1994; Berg and Scherg 1994; Pascual-Marqui et al. 1994). *NicoletBiomedicalInc., Madison, WI, U.S.A. ^Department of Electricaland Computer Engineering, University of Wisconsin - Madison, Madison, WI, U.S.A. #Department of Clinical Neurophysiology, University Hospital Utrecht, Utrecht, The Netherlands. Accepted for publication: April 16, 1996. The authors wish to thank Dr. J. M. Guerit, Dr. R.E. Lasky and G. Rook for their valuable suggestions. We thank Dr. K. Hecox for his support and R. Birrenkot for preparing the figures. Correspondence and reprint requests should be addressed to Dr. Wire van Drongelen, Nicolet Biomedical Inc., 5225 Verona Road, P.O. Box 44451,Madison, W153744-4451,USA. Copyright 9 1996Human SciencesPress, Inc.
All studies are based on a m o d e l that describes h o w an electromagnetically active area in the brain generates a field that can be detected b y electrical or magnetic sensors. Theoretically, this field can be p r e d i c t e d f r o m knowledge of the source locations, g e o m e t r y a n d electromagnetic properties of the source and m e d i u m . The problem of predicting the surface field from the source and m e d i u m is t e r m e d the f o r w a r d problem. O n the basis of the geometry of the head, skull and brain, a realistic solution can be obtained using a finite element, or b o u n d ary element, approach (Meijs 1988; H a m a l a i n e n and Sarvas 1989; Yan et al. 1991). Because these m e t h o d s are computationally intense, m a n y authors use a simplified model of the head, a sphere with three or four shells (Kavanagh et al. 1987; de M u n c k 1989; Salu et al. 1990). Some authors use a homogeneous sphere with a weighted dipole position and dipole value as an approximation of the multisphere m o d e l (Scherg and v o n Cram o n 1985; A r y et al. 1981; Berg and Scherg 1994). The localization p r o b l e m is the inverse of the forw a r d solution: the field is sampled and the location of the source needs to be determined. For the inverse problem, a unique relation b e t w e e n the r e c o r d e d field and source does not exist. The m e a s u r e d data can correspond to different source distributions. Additional criteria must be specified to determine a unique source distribution. Several approaches have been a d o p t e d (Scherg 1992; Z h a n g and Jewett 1994). Each has its particular constraints: e.g., the n u m b e r of active sources is guessed, or a fixed location and orientation of the dipoles over time is assumed
40
(Scherg 1992). A central issue in most localization algorithms is the estimation of the number of active sources at any given point in time. The present study evaluates a method for localizing electrical activity that is based on principles of spatial filtering. Data collected at multiple surface sites is combined in such a w a y as to pass source activity originating at a given location in the brain while attenuating activity originating from other locations. The covariance matrix of the measured data is used to design the spatial filter. The power in the filter output as a function of filter passband location provides an estimate of source activity power as a function of location within the brain. This approach does not require any assumptions about the number of active sources. Indeed, under the appropriate conditions, both the number and locations of sources may be determined from the map of estimated power as a function of location. The efficacy of this approach is studied using both simulated and measured data. Simulations were used to examine the relationship between the resolving power of the method, the number of recording sites, the signal-to-noise ratio and the depth. The measure d electrophysiological data were selected to test the method's ability to detect activity from different locations. The method we describe is an adaptation of techniques used in RADAR and SONAR applications (Van Veen and Buckley 1988). Theoretical considerations in applying the method to neurophysiology were published previously (Van Veen et al. 1992). Applications of spatial filtering to the localization problem have also been discussed by Mosher et al. (1992) and Spencer et al. (1992). A technique for estimating variance in source activity (Dale and Sereno 1993) employs similar concepts to those used here.
Methods and Material Spatial Filter In the following, the method is referred to as linearly constrained minimum variance (LCMV) filtering. Boldface lower and upper case symbols denote vectors and matrices respectively. The method is based on the concept of spatial filtering that parallels the more familiar concept of temporal filtering. In the latter filter one discriminates signals based on their temporal frequency content. Here we discriminate signals based on their spatial origin. A 'narrowband' spatial filter passes signals originating from a small 'passband' volume while attenuating those originating from other locations. This spatial filtering is based on source and head models that relate the underlying neural activity to the electromagnetic field at the surface. The LCMV method can be applied to
van Drongelen et at.
both EEG and MEG, however, in the text we will describe the model for electrical sources. Neural activity is modeled as a current dipole or sum of current dipoles. Let x be an N by 1 vector composed of the potentials measured at the N electrode sites The potential due to a single dipole with location repre sented by vector q is expressed as: x = H(q).m(q). Here the elements of the 3 by I vector m(q) are the X, Y and Z components of the dipole moment and the columns of the N by 3 matrix H(q) represent solutions to the forward problem. That is, the first column of H(q) is the potential at the electrodes due to a dipole source at location q having unity moment in the X direction and zero moment in Y and Z directions. Similarly, the second and third columns represent the potential due to sources with unity moment in Y and Z directions respectively. The matrix H will be referred to as the transfer function in the following. N o w suppose x is composed of the potentials due to L active dipole sources at locations Cli,i=1,2 ..... L and noise n. The medium is linear so the potential at the scalp is the superposition of potentials from the active neurons. Hence, L
x = E H(qi )"m(qi ) + n i=l
(1)
For very large L, the sum may be represented as an integral over the volume containing possible sources. Note that x does not contain any temporal information since it is obtained by sampling the electrodes at a single time instant. Since neural activity has a probabilistic component, we model the dipole moment as a random quantity and describe its behavior in terms of mean and covariance C(qi). We assume the noise is zero mean with covariance matrix Q and the moments associated with different dipoles are uncorrelated. The covariance matrix of measured potentials C(x) is then expressed as:
L C(x) = E H(qi )" c(qi )"HT(qi ) + Q
(2)
i=1
Where HT(qi) is the transpose of H(ql). The goal of the LCMV method is to design a bank of spatial filters where each filter passes signals originating from a specific location in the brain while attenuating signals from other locations. Define the spatial filter for a narrowband volume element Q0 centered on location q0 as the N by 3 matrix W(q0) and let the three component filter output y be the inner product of W(q0) and x: y = WT(qo).x
(3)
Spatial Filtering in Source Localization
41
An ideal narrowband filter satisfies:
f Iifq~Q0
W T ( q o ) ' H ( q ) = ~ ~ ~ ,, c~ [o ~ ~l = "~o and q ~ B
(4)
where B represents the volume of the brain. If (4) were satisfied, then in the absence of noise, the filter output is the dipole moment at the location of interest. As in temporal filtering, it is generally impossible to have complete attenuation in the stop band. Given this limitation, the question is h o w to design a tilter that in some sense is optimal. The LCMV approach offers a guiding principle for such a design. The idea is to find a W(q0) that minimizes the variance at the filter output while satisfying the linear response constraint based on (4), WW(q0).H(q0) = I. Hence the name linearly constrained minimum variance. The constraint insures that signals of interest are passed; while in minimization of variance, one minimizes the contributions of noise and energy outside the region of interest (q ~ Q0). We can summarize this constrained minimization problem as: min tr{C(y)} subjecttoWT(q0).H(q0)=I
(5)
W(q0) Where C(y) = WT(q0).C(x).W(q0) and tr{A} denotes the trace of matrix A. Using the method of Lagrange multipliers, the solution of (5) yields: W(qo) = (HT(qo).C'I(x).H(qo))q.HT(qo).C'I(x)
(6)
Equation (6) describes the LCMV spatial filter W as a function of the transfer function H and data covariance matrix C. Assuming the spatial filter passes signals originating at location q0 and attenuates signals from other locations, then substitution of (6) into the equation for the tilter output, equation (3), gives an estimate of the dipole moment at location q0. Under this assumption, the variance of the filter output tr{C(y)}, represents an estimate of the neural activity power associated with location q0Substituting (6) for W(q0) in the expression for C(y) gives an expression for the estimated variance as a function of location: Varest(qo) = tr{(HT(qO).C'l(x).H(qO)) q}
(7)
The subscript est denotes that this is an estimate of the neural activity power. By depicting the estimated variance as a function of location in the brain, we obtain an estimate of the neural power as a function of location. The variance estimate given in (7) represents neural
activity associated with sources and the effect of noise present in the data. It can be shown that noise introduces a spatially non-uniform component to the variance estimate with the actual noise contribution for a particular location depending on the transfer function matrix H. To prevent this non-uniform noise contribution from being mistaken for or masking actual neural activity, we form an activity index by normalizing the variance estimate (7) by an estimated noise activity level. That is, if the noise has covariance matrix Q, we divide (7) b y the estimated noise variance, tr{(HT.Q-1.H) q} to obtain the index. This activity index represents the ratio of the total activity to the normalized noise activity. Large values indicate regions of strong neural activity. The noise covariance matrix can be estimated from data that is collected in the absence of active sources. Alternatively we found that simply assuming Q = I often yields good results. This corresponds to assuming the noise in different electrode channels is uncorrelated and of equal variance. The absolute noise level is unnecessary since it simply scales the entire activity index by the same constant. In the simulations and experiments that follow, we depict the activity index as a function of location assuming Q = I. Simulation
The simulated data is based on a three-shell spherical model and calculated with the algorithm described by Salu et al. (1990). We assumed a headmodeI of 8.25 cm radius. The skull and scalp were 0.3 cm and 0.7 cm thick respectively. The ratio of conductivity between skull and soft tissue was set to 0.0125. Dipole moments were set to unity in the X, Y and Z directions. The X, Y and Z coordinates are front-back (front positive), left-right (left positive) and up-down (up positive) respectively. The model was used to produce both the transfer function H and to generate the potential values at the 'recording sites'. Transfer functions were calculated on a spatial grid having a resolution of 0.25 cm in the X and Y directions and a 1 cm resolution in the Z direction. Because the center value (0,0,0) generates computation problems calculations started at 0.125 for the X and Y directions. Two pairs of dipole sources with different depth were used to determine the resolving power of the method. One pair of sources was placed superficially, the other pair of dipoles close to the center of the sphere. The X, Y, Z coordinates were 1.375, 0.125,1.0 and -1.625, 0.125, 1.0 for the deep sources; for the sources close to the surface we used 1.375, 6.375, 1.0 and -1.625, 6375, 1.0. We also simulated 20,32 and 64 electrode arrays. The 'electrode' locations for which the potential is calculated are: 20 positions (Fpl Fp2 F7 F3 Fz F4 F8 T3 C3 Cz C4 T4 T5 P3 Pz P4 T6 O 1 0 z 02); 32 positions (Fp1 Fp2 AF3 AF4 F7 F3 Fz F4 F8 FC5 FC1 FC2 FC6 T3 C3 Cz C4 T4 CP5 CP1
42
CP2 CP6 T5 P3 Pz P4 T6 PO3 PO4 O 1 0 z 02); 64 positions (Fpl Fpz Fp2 AF7 AF3 AFz AF4 AF8 F7 F5 F3 F1 Fz F2 F4 F6 F8 FT7 FC5 FC3 FC1 FCz FC2 FC4 FC6 FT8 T9 T3 C5 C3 C1 Cz C2 C4 C6 T4 T10 TP7 CP5 CP3 CP1 CPz CP2 CP4 CP6 TP8 T5 P5 P3 P1 Pz P2 P4 P6 T6 PO7 PO3 POz PO4 PO8 O 1 0 z 0 2 Iz). Three different signal-to-noise conditions were studied by adding a different amount of white noise to the covariance matrix. The noise levels resulted in average signal-to-noise ratios of 95, 9.5 and 0.95. These noise levels indicate the ratio between the variance in the signal and the variance in the noise. The relation to the more conventional amplitude signal-tonoise ratio can be estimated as the square root of these values. The spatial filter was applied to the simulated 'recordings'. The results were depicted in 8 horizontal cross-sections through the spherical model. The lowest section intersects the center of the sphere, with the other slices separated b y I cm. Experimental d a t a
SEP recording. The SEP data was obtained by using two Nicolet Pathfinder I systems simultaneously. Each system received input from a 32 channel amplifier (SM2000); one system triggered the other to provide synchronization of the averages. In this manner 64 channels of data were recorded. The montage used is: Fpl Fpz Fp2 AF7 AF3 AFz AF4 AF8 F7 F5 F3 F1 Fz F2 F4 F6 F8 FT7 FC5 FC3 FC1 FCz FC2 FC4 FC6 FT8 T3 C5 C3 C1 Cz C2 C4 C6 T4 TP9 TP7 CP5 CP3 CP1 CPz CP2 CP4 CP6 TP8 TP10 T5 P5 P3 P1 Pz P2 P4 P6 T6 PO7 PO3 POz PO4 PO8 O 1 0 z 0 2 Iz. The average value of A1 and A2 is used as the reference. The filter had a passband between 5 and 250 Hz. The signals were digitized at a rate of 2560 samples per channel per second. Electrode locations were digitized with a 3Space Isotrak (Polhemus Navigation Sciences). The current stimulus was delivered as a pulse of 100 us at a rate of 7.2 stimuli per second and an intensity that elicited slight movement of the thenar muscle. The SEPs were performed on a healthy test subject. Three experiments were performed: left-, right- and bilateral median nerve stimulation. In each experiment 36 averages were acquired consisting of 200 trials each. The total timebase used was 100 ms with a 20 ms prestimulus interval. The covariance matrix was calculated using five points (an interval of 1.56 ms) around the N20 peak. An interval of the subaverages around the N20 peak was digitally filtered. The passband of the digital filter was set between 10 and 100 Hz. The average value of the prestimulus period was subtracted from the corresponding recorded waveform to correct for offset. A least squares fit of the measured electrode positions to a spherical surface was p e r f o r m e d and a three shell spherical
van Drongelen et at.
model employed for computing the transfer functions. The radius of the sphere was based on measured head size, the skin and skull thickness was set to 0.3 cm and 0.7 cm respectively. The conductivities were identical to the ones used in the simulations.
Spike data The spike data were recorded in the laboratory of one of the authors (vH). The data of the 2 patients presented were selected from the laboratory archive on the basis of multiple active sites, availability of additional structural data and depth recordings. The mon, tages used in both cases were different and adapted to the location of the epileptiform activity. In one patient (LH) the activity was located occipitally, the other patient (RY) showed centrotemporal activity. The two 32 channel montages are as follows: LH (Fpl F2 F7 F3 Fz F4 F8 FC1 FC2 T3 C3 Cz C4 T4 TP9 CP5 CP1 CP2 CP6 TP10 T5 P3 Pz P4 T6 POz P O 7 PO8 O1 0 2 Iz A1); RY (Fpl Fp2 AFz F7 F3 Fz F4 F8 FT9 FC5 FC1 FC2 F C 6 T 3 C3 Cz C4 T4 TP9 CP5 CP1 CP2 CP6 T5 P3 Pz P4 T6 O t POz 0 2 A1). In both montages, A2 is used as reference~ The EEG was recorded with a Nicolet Pathfinder II connected to a 32 channel amplifier (SM2000) with a band-pass filter set between 0.5 and 70 Hz and a sampling rate of 250 samples per channel per second (van der Meij 1992). In both cases, the covariance of the spike signal was calculated from single time series including spike activities (figures 3 and 4). The parameters (ge~ ometry and conductivities) of the model were identical to the ones in the simulations. Because the activity in the patients could be located in the lower parts of the brain, we included analysis of additional cross-sections at 1, 2 and 3 cm below the center of the sphere. The coordinates of the electrode locations were estimated on the basis of an ideal placement on a sphere.
Results Simulations
All the simulated dipoles were correctly located by the LCMV algorithm in the horizontal slice at Z=I cm. For this reason we will only present the results for this cross-section. The fall-off of the peaks in the neural activity index in the Z direction is of the same magnitude as in the X and Y directions. The results for the different parameters are summarized in the 3-D mesh plots in figure 1. The X and Y axes represent the coordinate values in the slice, the Z value in figure I is the activity index plotted versus X and Y. A dipole source is associated with a peak in the 'landscape' of the activity index.
S p a t i a l Filtering in S o u r c e L o c a l i z a t i o n
43
xlO 4
•
10 ~ 3-
go.sj
~ 2.5--~
2-
9~ Ls. ,~ 1. 0.5O5
6
5
6 0
u (cm)
X (era)
Y (cm)
x 10 ~
X (cm)
Y (cm)
x 10 4
X (cm)
10 ~ 2.5F 2,
O-
5 ~'--...~ 5 u ~ ~ ' ~ : ' ~ 0 .5~ > ~ _ _ ~ 6 Y (cm)
x
_4
~
2
'
0
6
~
-2 X (crn]
0
Y (crn)
104
X (crn)
Y (cm)
•10~ •
0.8
0.04
,
~o. 5
6
Y (era)
5
-6
X (cm)
I
G x (cm)
x10~
H
1.5
~ i
~
~o.
~o. 5
6
Y (cm)
X (crn)
5
-
Y (crn)
6
X (cm)
Figure 1. Results of a simulation study. Spatial filtering is used for localization of electrical dipoles in a three-sphere model. To study resolving power of the calculations, a pair of superficial and a pair of deep dipoles simulate electrical activity in the brain. All four dipoles are located in a section 1 cm above the center of the sphere. In each of the graphs, the output of the filter, the activity index, is plotted versus the X, Y location in a cross-section at Z=I cm. Spatial resolution decreases with distance from the surface. Comparing the graphs one can observe the effects of the number of electrodes: 64 positions (A, B, C); 32 (D, E, F) a n d the 20 positions of the 10-20 system (G, H, D. The effects of the signal-to-noise ratio c a n be established by comparing A, D, G (average ratio = 95) with B, E, H (average ratio = 9.5) and C, F, I (average ratio = 0.95).
The set of parameters with maximal resolving power of the examples presented is the low noise, 64 channel simulation depicted in figure 1A. The two pairs of dipoles were located correctly: i.e., the maxima of the four peaks are located at the source locations. In this solution, the peaks of the deep pair of dipoles are smeared relative to the superficial ones and they merge in the center of the plot. The superficial pair of sources appear as separate and sharp peaks. When the noise level is increased for the 64 channel situation, the separation between the central peaks disappears (figure 1B and 1C). The two separate superficial dipoles are diffi-
cult to distinguish in the highest noise level (figure 1C). If the number of channels is decreased to 32, the solution for the low noise situation shows four peaks with the maxima at the correct locations (figure 1D). Increasing the noise causes the two deep peaks to merge (figure 1E). In the lowest signal-to-noise simulation one can only distinguish between deep and superficial activity (figure 1F). The last set of simulations was performed with the 10-20 montage. Although the peaks are smeared and merged, correct answers were obtained for the maxima in the low noise calculation (figure 1G). The solution in figure 1H distinguishes two regions of
44
van Drongelen et al.
Left
Frontal
Z
model, therefore the activity index is mapped into a sphere. Accordingly, we will refer to anatomical correlates of the locations in a global manner. A coarse correlation between our localized position and anatomy can be made if one refers to the electrode positions. In the simulations and the analyses of the spike data, the FpI, Fp2, and occipital electrodes are located slightly below the cross-section through the center of the sphere at Z=0 (0.3 cm below Z=0 in a sphere with a radius of 8 cm); Cz is located at X,Y,Z = 0,0,R with R representing the radius of the sphere. For the SEP data analysis, the measured positions deviated only slightly from the idealized ones used in the other localizations.
SEP Right Left
B
Frontal
\ Z
Right Figure 2. Localization of the N20 in the SEP upon bilateral median nerve stimulation. Maxima in the neural activity index were found 4-5 cm a b o v e the center of the sphere. The small half sphere diagrams represent a coronal section (the YZ plane) through the head model with an indication of the level of the contour plots of the activity index 'landscape' in A and B (4 and 5 cm above the center of the sphere respectively). The two sources, located left and right in the posterior half of the cross-sections are compatible with two active areas in the left and right sensory cortex.
activity. In the simulation with the worst conditions in the set of examples, there is only an indication for unilateral activity in the cross-section (figure 1I). Experimental d a t a The forward solutions were based on a three-shell
The results of bilateral stimulation of the median nerve from one of the subjects are shown in figure 2. Unilateral left and right stimulation was used to confirm the localization of the early cortical activity in the con: tralateral hemisphere. In both cases the N20 activity was localized as a single source in the correct hemisphere that corresponds well with the sensory cortex given the resolution of the spherical model. Two simultaneously active sources were expected upon bilateral stimulation. The main question was whether the resolving power of the algorithm was sufficient to distinguish these two sources. In the data obtained during bilateral stimulation we found a left and a right peak located at 4-5 cm above the center of the sphere. Figures 2A and 2B depict the activity index at 4 cm and 5 cm above the center of the sphere respectively. The peaks were found in the same locations as the individual peaks found in the data of the unilateral SEP study.
Occipital spikes In the subject identified as LH, an enlarged left ventricle was detected prenatally (figure 3A). LH showed a slight psychomotor retardation, a mild right sided hemiparesis and a h o m o n y m o u s hemianopsia. Seizures started at 1.5 years of age. Convulsions consisted of loss of consciousness and jerking of the right bodyhatf. Postictally there was a transient increase of hemiparesis. Seizures frequently progressed to a status epilepticus. The seizures became resistant to drug therapy and it was decided to perform a leftsided hemispherectomy at 6 years of age. This surgery resulted in a cessation of the seizures and a slight increase of the hemiparesis. A typical recording from electrode position POz is shown in figure 3B. We examined the localization of the two larger spikes indicated with I and 2 in figure 3B. The covariance matrix was based on 20 samples of each of the spikes. A contour plot of the activity index 'landscape! is shown in figure 3C. Two larger sources were detected in
Spatial Filtering in Source Localization
45
t2
I!/ I c
,
Frontal
Left
~ ~
Right
Occipital Figure 3, A summary of the results obtained from the epilepsy data from patient LH. An MRI is depicted in A; the enlarged lateral ventricle surrounded by a thin layer of cortex dominates the left hemisphere. Note that the left side of the brain is on the right side in the image. The recording in B is a typical epoch of the EEG at position POz. The two spikes labeled 1 and 2 were used to calculate the covariance matrix. Calibrations in B indicates 0.2 s and 50 uV. The active locations associated with this data are shown in the contour plot of the activity index in C, Two locations (arrows) are found in the left hemisphere: one source is located in the lateral part, the other is located close to the midline.
the section at 1 cm below the center of the sphere, the slice that is depicted. A comparison with the structural data in figure 3A indicates that these sources may represent activity in the medial and lateral parts of the thin layer of occipital cortex that surrounds the enlarged ventricle. Indeed, electrocorticography during surgery showed abundant spiking in the left occipital region, spreading towards temporal and parietal areas. A strip electrode positioned interhemispherically revealed continuous spiking at the medial surface of the occipital lobe.
Eye movement Many eye movement artifacts were found in the EEG
of patient RY. Often the movement artifacts coincided with spike activity, an example is shown in figure 4A. The Fpl and Fp2 electrodes show a large artifact, whereas the spike activity in T3(T7) is lower in amplitude. In several EEG epochs, we used combined eye artifact and spike activity to calculate a covariance matrix which included the activity of at least three sources: both eyes and at least one spike source. In all of the cases, we could distinguish the centro-temporally located spike from the eyes. In many cases the eyes were identified as a single source located at the median between the real positions. This is not an unexpected finding because eye movements produce highly correlated electrical signals and the method separates uncorrelated activity only. In some
46
van Drongelen et ol.
A
I
t I.,I I I,I
,x tcmj
X
' ~"'V
A (C/T1)
Figure 4, An example of spatial analysis of a combined eye movement artifact and spike activity (arrow in A), in A, the EEG ofT3 (upper trace), Fpl (middle) and Fp2 (lower) is shown, Calibration lines indicate 0,2 s and 50 uV. The EEG interval between the vertical lines in A w a s used to calculate the covariance matrix, Figures B o n d C are 3D mesh plots of the activity index plotted versus X and Y in cross-sections at 3 and 4 cm above the center of the sphere. The X coordinate is front-back (front-positive) and Y is left-right (left-positive). The 'landscape' in C shows a maximum associated with the spike activity (s): no activity is seen at the locations of the eyes. Some of the spike-actMty is also present in B, a plot of the activity index in the section 1 cm lower. The two lower peaks (e) left and right at the frontal rim in plot B are associated with the eye movement signal.
cases both eyes were located in a 'real' left and a right source, an example can be seen in figure 4.13igures 4B and 4C show the activity index obtained from the data shown in figure 4A. The part of the signal between the vertical lines in figure _4Aincludes a large amplitude artifact and a small spike. In the 3D mesh plot of figure 413, the maximum activity index represents the spike; the two smaller maxima at the frontal rim of the plot correspond with the eye locations. Spatial analysis of the covariance matrix of this epoch shows a high activity in the lower cross-sections (maximal at 2 cm below the center of the sphere) probably corresponding to the eye movement artifact and a peak of activity at 3-4 cm above the certter of the sphere. In spite of the large amplitude difference, we find the largest activity index values at the centro-
temporal location in the higher cross-sections at 3 and 4 cm above the center. The positions of both eyes show slight activity at 3 cm (figure 4B); the activity of these positions disappeared at 4 cm (figure 4C) where the spike location showed the maximum of activity,
Discussion It is well known that there is no unique solution to the electromagnetic inverse problem. For this reason all proposed models for localization of neural activity rely on one or several constraints. Ideally one would hope to find a set of constraints that cart be applied in all types of measurements, Alternativel}5 one could try to establish constraints that apply in a specific situation. In terms of
Spatial Filtering in Source Localization
the output of the source localization algorithms, two general approaches are currently used: equivalent dipole modeling (e.g., de Munck 1989; Scherg 1990) and tomographic reconstruction methods (e.g., Ioannides et al. 1993; Pascual-Marqui et al. 1994). In dipole modeling the number of activities, the geometry and electromagnetic parameters have a major impact on the solution. Most of the tomographic solutions are attractive because they do not require a priori estimation of the number of active sites. These solutions all resolve the current density distribution as a function of location. In the LORETA (Low Resolution Brain Electromagnetic Tomography) algorithm (Pascual-Marqui et al. 1994) the smoothest of all solutions is chosen as the output. In minimum norm algorithms, such as described by Ioannides et al. (1993), the sum of the square of all neural currents is minimized. The assumption in LORETA can be realistic in long latency event related potentials where larger areas in the brain may be active. The reahty of the minimum norm has not been tested, but tends to generate dispersed and superficial solutions (Gorodnitsky et al. 1995). Application of LORETA and minimum norm methods in recordings of focal activities will fail to generate realistic solutions. In all algorithms, including the one presented in this paper, the solution is sensitive to the model parameters used for the forward solution. This is good because one can apply the method for different recording-modalities and different case individual geometries. It is, however, also in, practical because a realistic solution requires computational intense algorithms. Application of high resolution EEG to obtain a map of the cortical activity may be a practical compromise to have a first estimate of the location and number of active cortical areas (Nunez et al. 1994). The LCMV method requires inversion of the covariance matrix. For this reason, the number of data vectors (one vector being represented by x in equation (1)) must be larger than the number of electrodes. One may estimate the covariance matrix from single points in repetitive observations, s u b s e q u e n t samples in a single recording, or a combination of both approaches. Several other current dipole localization methods do not have this requirement; these algorithms can generate solutions on the basis of a single data vector. These solutions may not be realistic because one must guess the number of active sources at the time instant being studied. Authors that apply linear decomposition require a sequence of sampled points and must assume stationary dipoles, an assumption that can be avoided using the spatial filtering approach described in this study. A method described as multiple signal classification (MUSIC) resembles the LCMV filter approach in that it generates a type of spatial activity index (Mosher et al.
47
1992; Spencer et al. 1992). However, MUSIC requires specification of the dimension of the signal subspace, which is analogous to requiring k n o w l e d g e of the number of sources. Dale and Sereno (1993) estimated dipole strength as a function of location over the entire interior of the brain. This is accomplished by minimizing the expected error between the estimated and corrected solution. This method requires knowledge of the dipole strength vector covariance. Their estimation of the dipole strength vector covariance is similar to the one described in this paper. The simulation results indicate excellent performance of the LCMV algorithm when the signal-to-noise is high (figure 1A, D, G). The number of electrodes has a similar effect: more electrodes are associated with a higher resolving power. Finally, resolving power decreases with depth. The simulation study indicates that the LCMV method may be expected to perform well with cortical sources that are situated superficially. Deeper activity can be detected, but individual sources may be lost because of smearing (e.g., figure 1B, C, E, F, H, I). The forward solution in the simulations can be matched with the source locations because the same model generates the transfer function and the simulated recorded data. In reality a three-shell model only approximates the real field and can generate considerable error (Roth et al. 1993). Furthermore, real activity may exhibit considerable correlation between active sites and the noise characteristics may differ from that assumed in the simulations. For these reasons, the performance of the algorithm applied to real recorded data can not be expected to exceed that observed in the simulations. From a simulation study we performed, we know that violation of the assumptions reduces the resolving power to some extend. Because it is beyond the scope of this paper, quantitative effects of noise distribution and correlated activity have not been addressed. Application of the algorithm in noisy data recorded with only 20 electrodes can, at most, indicate areas of activity (figure 1H, I). When the algorithm cannot resolve the locations of activity, a smeared peak of neural power is obtained. In many situations electrophysiological data cannot be used to verify localization methods. The precise location of the measured activity is simply not known and verification is impossible. In some cases, outcome can be falsified: sources in a ventricle, the longitudinal fissura or even outside the head are examples. In our case we chose the N20 of the median SEP because of the consensus as to its source in the contralateral parietal cortex (Mauguiere et al. 1983; Allison et al. 1989; Franssen 1990). On the basis of the results depicted in figure 1, we can compare the SEP recording with the simulations. Because the amplitude signal-to-noise varies from I to 10, we expect to be in a situation in between the simulations shown in figures
48
1A and lB. With the N20 being a source with a relatively superficial location, one m a y expect a correct location and sufficient resolving power to distinguish between the right and left peak in the bilateral stimulation experiment. The SEP test (figure 2) shows that the algorithm is indeed able to detect the active areas and separates the sources in a bilateral stimulation experiment. Because of the approximation of the head with the three-shell model, exact anatomical correlates cannot be obtained. However, the locations and s y m m e t r y of the sources left and right are in good agreement with the expectation of two active sources in the sensory cortex of both hemispheres. We used 32 scalp electrodes for the analysis of the spikes. Because of the high signal-to-noise in this type of signal, we assume that this situation resembles the simulation in figure 1D or 1E where superficial sources can be seen separately. The spike data was selected because the MRI of one patient showed a thin layer of cortex around an enlarged ventricle in the left hemisphere. From depth recording and clinical follow up, it is apparent that the epileptiform activity originated in this sheet as confirmed by the result in figure3C. We found the analysis of this situation challenging because simultaneous activity in the thin curved cortical layer would make a mis-localization of an equivalent source in the cerebrospinal fluid of the ventricle highly probable: i.e., the vector s u m of the two dipoles in figure 3C is situated in between both the detected locations in the enlarged ventricle. The other epilepsy case showed a high occurrence of combined eye m o v e m e n t and spike discharges (figure 4). The LCMV filter performance during a relatively high amplitude eye-movement artifact indicates the potential of this approach to discriminate spatially distinct signals. The result shown in figure 4C indicates that the spatial filter attenuates strong signals coming from the eye, while the location of the relatively weak spike activity displays a high activity index.
References Allison, T., McCarthy, G., Wood, C.C., Darcey, T.M., Spencer, D.D. and Williamson, P.D. Human cortical potentials evoked by stimulation of the median nerve. I. Cytoarchitectonic areas generating short-latency activity. J. Neurophysiol., 1989, 62:694-710. Amir, A. Uniqueness of the generators of brain evoked potential maps. IEEE Trans. Biomed. Eng., 1994, 41:1-11. Ary, J.P., Klein, S.A. and Fender, D.H. Location of sources of evoked scalp potentials: correction for skull and scalp thickness. IEEE Trans. Biomed. Eng., 1981, BME-28: 447452. Berg, P. and Scherg, M. A fast method for forward computation of multiple-shell spherical head models. Electroenceph. Clin. Neurophysiol., 1994, 90:58-64. Buchsbaum, M.S., Hazlett, E., Sicotte, N., Ball, R. and Johnson,
van Drongelen et ai.
S. Geometric and scaling issues in topographic electroencephalography. In: F.H. Duffy (Ed.), Topographic mapping of brain electrical activity, Butterworths, Boston, 1986: 325-337. Chatrian, G.E., Lettich, E. and Nelson, P.L. Ten percent electrode system for topographic studies of spontaneous and evoked EEG activities. Am. J. EEG Technol., 1985: 83-92. Dale, A.M. and Sereno, M.I. Improved localization of cortical activity by combined EEG ancl MEG with MRI cortical surface reconstruction: a linear approach. Journal of Cognitive Neuroscience, 1993, 5:162-176. Ebersole, J.S. EEG dipole modeling in complex partial epilepsy. Brain Topography, 1991, 4:113-123. Fender, D.H. Source localization of brain electrical activity, in: A.S. Gevins and A. Remonds (Eds.), Methods of Analysis of Brain Electrical and Magnetic Signals. EEG Handbook, revised series Vol. 1, Elsevier, Amsterdam, 1987: 355-403. Franssen, H. Relation between topography and generators of somatosensory evoked potential (SEP) to median nerve stimulation. Neuro-Orthopedics, 1990, 9: 77-88. Gorodnitsky, I.F., George, J.S. and Rao, B.D. Neuromagnetic source imaging with FOCUSS: a recursive weighted minimum norm algorithm. Electroenceph. Clin. Neurophysiol., 1995, 95:231-251. Hamalainen, M.S. and Sarvas, J. Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data. IEEE Trans. Biomed. Eng., 1989, 36:165-171. Ioarmides, A.A., Singh, K.D., Hasson, R., Baumann, S.B., Rogers, R.L., Guinto, F.C. and Papanicolaou, A.C. Comparison of single current dipole and magnetic field tomography analysis of the cortical response to auditory stimuli. Brain Topography, 1993, 6:27-34. Kavanagh, R.N., Darcey, T.M., Lehman, D. and Fender, D.H. Evaluation of methods for three-dimensional localization of electrical sources in the human brain. IEEE Trans. Biomed. Eng., 1987, BME-24:421-429. Mauguiere, F., Desmedt, J.E. and Courjon, J. Astereognosis and dissociated loss of frontal or parietal components of somatosensory evoked potentials in hemispheric lesions. Brain, 1983, 106:271-311. Meij, van der W. Rolandic epilepsy. Clinical significance of quantitative spatial and temporal EEG analysis of rolandic spikes. Thesis, Department of Clinical Neurophysiology, University Hospital Utrecht, Utrecht, 1992. Meijs, J.W.H. The influence of head geometries on electro- and magnetoencephalograms. Thesis, Biomedical Engineering Division, University of Twente, Enschede, 1988. Mosher, J.C., Lewis, P.S. and Leahy, R.M. Multiple dipole modeling and localization from spatio-temporal MEG data. IEEE Trans. Biomed. Eng., 1992, BME-39:541-557. Munck, de J.C. A mathematical and physical interpretation of the electromagnetic field of the brain. Thesis, Ophtalmic Research Institute, University of Amsterdam, Amsterdam, 1989. Munck, de J.C. and Peters, M.J. A fast method to compute the potential in the multisphere model. IEEE Trans. Biomed. Eng., 1993, 40:1166-1174. Nunez, P.L. Methods to estimate spatial properties of dynamic cortical source activity. In: G. Pfurzscheller and F. Lopes da
Spatial Filtering in Source Localization
Silva (Eds.), Functional Brain Imaging, Hans Huber, Toronto, 1988: 3-10. Nunez, P.L., Silberstein, R.B., Cadusch, P.J., Wijesinghe, R.S., Westdorp, A.F. and Srinivasan, R. A theoretical and experimental study of high resolution EEG based on surface Laplacians and cortical imaging. Electroenceph. Clin. Neurophysiol., 1994, 90:40-57. Pascual-Marqui, R.D., Michel, C.M. and Lehman, D. Low resolution electromagnetic tomography: a new method for localizing electrical activity in the brain. Int. J. of Psychophysiology, 1994, 18:49-65. Roth, B.J., Balish, M., Gorbach, A. and Sato, S. How well does a three-sphere model predict positions of dipoles in a realistically shaped head? Electroenceph. Clin. Neurophysiol., 1993, 87:175-184. Saiu, Y., Cohen, L.G., Rose, D., Sato, S. Kufta, C. and Hallett, M. An improved method for localizing electric brain dipoles. IEEE Trans. Biomed. Eng., 1990, 37:699-705. Scherg, M. Fundamentals of dipole source potential analysis. In: F. Grandori, M. Hoke and G.L. Roman (Eds.), Auditory evoked magnetic fields and electric potentials. Advances in Audiology, Vol. 6, Karger, Basel, 1990:40-69. Scherg, M. Functional imaging and localization of electromagnetic brain activity. Brain Topography, 1992, 5:103-111. Scherg, M. From EEG source localization to source imaging. Acta Neurol. Scan&, 1994, Supplement 152:29-30. Scherg, M. and von Cramon, D. Two bilateral sources of the late
49
AEP as identified by a spatio-temporal dipole model. Electroenceph. Clin. Neurophysiol., 1985, 62:32-44. Scherg, M. and Ebersole, J.S. Brain source imaging of focal and multifocal epileptiform EEG activity. Neurophysiol. Clin., 1994, 24:51-60. Spencer, M.E., Leahy, R.M., Mosher, J.C. and Lewis, P.S. Adaptive filters for monitoring localized brain activity from surface potential time series. Twentisixth Annual Asilomar Conference on Signals, Systems and Computers Conference Record, November 1992:156-161. Van Veen, B.D. and Buckley, K.M. Beamforming: a versatile approach to spatial filtering. IEEE ASSP Magazine, 1988, 5:4-24. Van Veen, B.D., Joseph, J. and Hecox, K. Localization of intracerebral sources of electrical activity via linearty constrained minimum variance spatial filtering. IEEE Workshop Statistical Signal and Array Processing proceedings, October 1992: 526-529. Wong, P.K.H. and Weinberg, H. Source estimation of scalp EEG focus. In: G. Pfurtscheller and Lopes da Silva (Eds.), Functional Brain Imaging, Hans Huber, Toronto, 1988: 89-95. Yan, Y., Nunez, P.L. and Hart, R.T. Fipdte-element model of the human head: scalp potentials due to dipole sources. Med. & Biol. Eng. Comput., 1991, 29:475-481. Zhang, Z. and Jewet, D.L. DSL and MUSIC under model misspecification and noise-conditions. Brain Topography, 1994, 7:151-161.
