WO2006012913A1 - Method and apparatus to localize, image and analyze discrete sources and their activity patterns - Google Patents

Abstract

A method of determining a finite number of discrete sources, in particular electrical or magnetic sources within a 3-dimensional volume as eg the human or animal brain, comprises the steps of: measuring a physical quantity using a number of n detectors having a known position, dividing the 3-dimensional volume into a plurality of discrete compartments, selecting a subset of m compartments, determining estimated activity values (Q-values) of the m compartments based on the measurement data, storing and updating the estimated activity values, iteratively repeating steps c), d) and e) with new subsets of compartments and, determining a number of compartments as active sources based on the stored and updated estimated activity values.

Description

Method and apparatus to localize, image and analyze discrete sources and their activity patterns in a 3D volume

Field of the invention

The present invention relates generally to a method of estimating source activity within a 3-dimensional (3D) physical volume from the data measured by a number of detectors on a 2-dimensional (2D) surface over or within the 3-D volume. In particular, the invention may be used for the localization of electrical and magnetic brain activity from the data obtained by electroencephalography (EEG) and magnetoencephalography (MEG), respectively. These techniques measure the electrical potential distribution (EEG) or the magnetic flux distribution (MEG) on a 2D-surface (the scalp) around a 3D-volume (the brain).

Background of the invention

Physicians and researches often wish to identify regions of electrically active cortex or other tissue in order to identify a source of illness or to map and monitor brain activity. While known monitoring equipment is capable of determining that electrical or magnetic activity has occurred, the location, strength and orientation of the source activity over time must be calculated. This calculation requires a model which relates the neural activity at a certain position within the volume to the distribution of electric potential or magnetic field on the surface. The most commonly used head-model for this so-called "forward" problem is a 3-shell spherical model.

The present invention is directed to solving the "inverse" problem, i.e. the calculation of the intra-cranial source activity responsible for an observed potential or magnetic field at the scalp surface.

One of the problems of estimating electrical and magnetic brain activity from such observations is the fact that it is an underdetermined problem: The brain consists of billions of neurons carrying electrical activity, whereas EEG and MEG are recorded at relatively few (21 to 306) locations. The data obtained from 21 recording channels is only sufficient to identify 7 independent sources, since each source is represented by a 3D current vector. A higher spatial resolution of source activity may only be obtained if some constraints and assumptions are made about the source activity.

A number of methods to solve the inverse problem of EEG and MEG have been proposed in the past, all using different assumptions. Most researchers assume minimum norm criteria to minimize the total power of the source currents in order to be able to estimate a very large number of small current dipoles throughout the brain or along the cortical surface of the brain. Others request maximum smoothness for the current source density in the brain. Both boundary conditions violate basic knowledge on the physiological activities of the brain (Scherg and von Cramon 1985, 1986; Scherg 1992; Scherg and Berg 1996, incorporated herein by reference).

An overview of various single-time-point and spatio-temporal methods may be found in Michael Scherg: "Functional imaging and localization of electromagnetic brain activity", brain topography, vol. 5, No. 2, 1992, also incorporated herein by reference. This paper also introduces Regional Source Imaging (RSI), a method based on a model of multiple dipoles with specific local constraints. A regional source consists of three collocated orthogonal dipoles. Successively one, two, three and more regional sources are used to "scan" the brain, the volume of which can be divided up into a set of discrete voxels. At the end of the scanning process, those voxels are selected for which each regional source images a maximal amount of variance while having the least amount of covariance with the other regional sources. If the number of regional sources exceeds the number of underlying macroscopic areas of activity, one or more regional sources begin to show source waveforms of low amplitude. This indicates that there is little contribution from these brain regions while there is localized contribution from the active regional sources. The inactive sources act like probe sources which check the surrounding brain region for the presence of source currents.

Another approach to the inverse problem is disclosed in US 5,263,488. This method uses a bank of spatial filters. Each filter essentially passes signal energy from one specific location within the brain while rejecting signal energy from other locations within the brain. This is achieved by a linearly constrained minimum variance criterion. The outputs of the filters can be used to estimate the signal power or dipole moment at each location, and this information can be displayed on a display device to provide a map of source activity within the brain. The spatial filters do not require prior knowledge of the number of sources, and the number of discrete sources can be identified as well as the location, power, and dipole moment of the sources. A disadvantage of the method described in US 5,263,488 is that only one voxel at a time can be processed. Further, the method of this document assumes that the activity of the discrete sources is described by orthogonal functions. However, since the neuronal activity in different regions of the brain is often correlated in space and time, this assumption may not hold.

Hence, there is a continued need to provide methods for estimating source activity within a 3-dimensional physical volume which do not require assumptions about the underlying source activity.

Summary of the invention

It is therefore an object of the invention to provide an improved method of estimating source activity within a 3-dimensional physical volume from the data measured by n detectors on a 2-dimensional surface over or within the volume, which is fast, may be completely automated and does not require assumptions on the orthogonality of the source signals, on the minimum norm of the source signals, or on maximum smoothness of source current density.

It is a further object of the invention to provide a method of imaging a 3-dimensional volume based on a finite number of discrete sources.

It is a still further object of the invention to provide a method of analyzing signals emitted by discrete sources in a 3-dimensional allowing an improved detection of pathological patterns.

The first object is achieved by a method of determining a finite number of discrete sources within a 3-dimensional volume, comprising the steps of measuring a physical quantity using a number of n detectors having a known position, dividing the 3- dimensional volume into a plurality of discrete compartments, selecting a subset of m compartments, determining an estimated activity value (Q-value) of the m compartments based on the measurement data, storing and updating the estimated activity values, iteratively repeating the last three steps with new subsets of compartments, and determining a number of compartments as active sources based on the stored and updated estimated activity values.

The second object is achieved by a method of imaging a 3-dimensional volume utilizing a finite plurality of discrete sources, comprising the steps of dividing the 3-dimensional volume into a finite plurality of compartments, providing locations and/or orientations of a defined number of active sources within the 3-dimensional volume, measuring a physical quantity at a plurality of n detectors, determining the activities of the active sources based on the measurement results, and creating a 3-dimensional image of the 3-dimensional volume by calculating the activity for each compartment using a probe source in addition to the active sources.

The third object is achieved by a method of time and frequency dependent pattern analysis of the activities of a plurality of discrete sources within a 3-dimensional volume, comprising the steps of dividing the 3-dimensional volume into a finite plurality of compartments, providing locations and/or orientations of a defined number of active sources within the 3-dimensional volume, time-resolved measuring a physical quantity at a plurality of n detectors, determining activities of the active sources based on the measurement results, deriving parameter vector representations of the three parameters time of occurence, spectral content, and duration of the activity of each of the active sources, performing a spectral-temporal filtering operation by comparing the parameter vector representations of the measured source activities with a previously stored two- dimensional set of basic signal patterns, or with a previously determined two-dimensional set of signal patterns derived from the complete measurement over time using spectral- temporal decomposition techniques on the source signals, respectively, and performing an iterative optimization procedure to determine as analysis result a set of discrete activity events for each active source defined by location, orientation, time of occurence, duration and periodicity or signal frequency, respectively.

Brief Description of the Drawings

Further objects, features and advantages of the invention will be apparent from the following detailed description when taken in conjunction with the accompanying drawings.

Fig. 1 shows illustrative pictures and a schematic block diagram depicting the basic principles of the present invention.

Fig. 2 is a schematic block diagram of the Discrete Source Scan Procedure according to the present invention.

Fig. 3 is a schematic block diagram of the Multiple Source Imaging Procedure according to the present invention.

Fig. 4 is a schematic block diagram of the Time and Pattern Definition Procedure according to the present invention. Fig. 5a is a schematic representation of the underlying physical model of a volume divided into a mesh of compartments and a plurality of detectors attached to an outer surface.

Fig. 5b shows the mathematical relationship between detector signals and source activity.

Fig. 6a is a schematic representation of the interactions between different sources of activity during external stimulation.

Fig. 6b shows examples of source activity signals.

Fig. 6c is a schematic block diagram illustrating the interactions between different sources of activity during external stimulation.

Fig. 7a is a schematic representation of the correlation between different sources of activity without external stimulation.

Fig. 7b shows examples of source activity signals.

Fig. 7c is a schematic block diagram illustrating the correlation between different sources of activity without external stimulation.

Fig. 8 is a more detailed schematic block diagram of the Time and Pattern Definition Procedure according to the present invention.

Figs. 9 and 10 are schematic block diagrams of different alternative implementations of the method according to the present invention.

Fig. 11 shows an exemplary hardware implementation of the present invention.

Detailed Description

The invention has been made for the processing of human EEG and MEG data, but it is not limited to electric potentials or magnetic fields recorded from the brain, heart or body. In fact, the invention can similarly be used to estimate the strength, localization and activity pattern of sources within a 3-dimensional (3D) volume from surface or depth measurements at n detectors in any physical system for which the propagation function of local source activities to the detectors is known, either by physical law or experimental observation. For example, the invention can be used to localize sources of heat, sound waves, magnetic disturbances, or centers of geological activity such as earthquakes from data taken at several points on the surface of the earth or another physical body. The main goal of the invention is to provide a method for automated and rapid localization and imaging of electrical activities in the human brain or body. Thus, the method and apparatus can be used to diagnose pathological abnormalities or monitor and reveal specific physiological functions.

Fig. 1 illustrates the principle of the invention using an example of an averaged EEG spike of a child with benign epilepsy. Electric potentials were recorded at 32 electrodes distributed over the scalp. Repeated spike events were recorded on the 2-dimensional (2D) scalp surface and averaged over an epoch of 500 ms. The discrete source scan procedure revealed 3 discrete centers, the 1^st & 2^nd in the left hemisphere in the vicinity of the central gyrus, and 3^rd at the homologous region in the right hemisphere. The topographical maps of propagated activity at the scalp surface due to the 3 sources, as given by the volume conductor model of the head, form the transfer matrix and are inverted to generate the inverse spatial filter (T^'1). Applying this spatial filter to the detector data matrix D yields the source activity signal matrix S. The 3 signals contained in S show spike patterns that overlap in time but propagate from source 1 to sources 2 and 3 as can be seen from the increasing time lag η of the spike peaks i=1..3. Using an additional probe source, the 3D functional image on the right can be created to reveal the source origins in a standardized (as depicted) or individual magnetic resonance image (MRI) of the brain.

The flow diagram in Fig. 1 highlights the major parts of the invention in gray. After preparation of the spatial data vectors, sources are found by the 'Discrete Source Scan Procedure'. The data vectors and the discrete sources are then used to create an image output by the 'Multiple Source Image Procedure'. The source image can be compared or fused with functional images obtained by other measurement modalities, for example, from functional MRI (fMRI) or positron emission tomography (PET). The detector data can also be processed in the time and frequency domains to reveal the specific times of occurrence and the patterns of the source activities by the Time & Pattern Definition Procedure'. Finally, the spatial, temporal and spectral parameters can be optimized in conjunction to define source activity events. The activity events can be determined, in principle, using a discrete mesh of compartments or voxels in a higher-dimensional space comprising the compartments of the 3D spatial volume plus a 3D pattern space spanned by the variables time of occurence, duration, frequency and phase, or by time of occurence, duration and periodicity of a set of basic patterns. Using additional parameters to describe the shape of the source activities, the 3D pattern space can be extended to higher dimensions. The first method of the invention, named 'Discrete Source Scan Procedure', may serve to calculate a 3-dimensional image of source activity over the whole 3D-volume on a finite mesh of compartments or voxels. The method assumes that the source activity in the brain may be represented by a finite number of discrete sources. In the brain, this assumption is generally valid, since neuronal activity occurs in the form of hundreds of thousands of neurons firing simultaneously in one circumscribed brain region (Scherg 1984). Thus, each active brain region can be modeled as one discrete equivalent source of electrical activity (Scherg 1984; Scherg and von Cramon 1985). The source model can be a single equivalent dipole at location r, with orientation angles 0, and ø, or a regional source with 3 orthogonal (or oblique) dipoles at center location r, that describes the source current vectors of any orientation in the surrounding brain region or compartment (Scherg and von Cramon 1986; Scherg and Berg 1996).

The main principle of the method according to claim 1 , as detailed below under 'Discrete Source Scan Procedure' and as shown in Fig. 2, is not to search only for the maxima of source activity in the 3D-volume - as is usually done - but rather to try out different subsets of voxels and determine a value of estimated source activity, called Q-value, for the whole set. Different methods of calculating the Q-value will be described below. The meaning of the Q-value is as follows: The Q-value depends on the locations of all voxels in the chosen subset and describes how much of the detected signal may be allocated to each voxel in the set. If the locations of the voxels of the set happen to coincide with the regions of activity, they will all exhibit a high Q-value. If some voxels of the set are near regions of activity, and the others are remote from regions of activity, the former will show a high Q-value and the latter will show a low Q-value. If all voxels are remote from regions of activity, they will all show a medium Q-value. Thus, the Q-value of a voxel is dependent on the position of the voxel and on the positions of the other voxels in the set. However, a voxel very near an active region will always display a relatively high Q-value and may thus be identified as an "active voxel".

According to one embodiment of the invention, the number m of voxels in the subset and the positions of the voxels in the subset are chosen randomly in each iteration. The Q- value may be computed very quickly so that a large number of sets of voxels may be tried within a short time. However, different strategies are defined in the dependent claims to speed up the procedure by using the information from previous iterations to choose selected voxels for the subsequent iterations additionally to the randomly placed voxels. One such strategy consists of excluding any voxels from further evaluation if their Q-value is below a predetermined or adaptive threshold criterion (KQ). Similarly, the voxels in the voxel set having the smallest Q-values, for example the lowest 5%, may be excluded from further evaluation. Since progressively more and more voxels are excluded, the procedure will get faster with each iteration. Preferably, it will stop a) when the number of remaining voxels, identified as the active voxels, is less than n. Alternatively, it will stop b) if the number of remaining voxels has not changed during the previous iterative steps, c) if a maximum number of iterations has been performed, or d), if a desired number of peaks with local maxima has been found (Fig. 2).

According to a further embodiment, it is possible to use some prior knowledge of the expected activity, i.e. the expected pattern of brain activity in response to a certain stimulus, to choose the positions of the voxels in the set. This knowledge may come, for example, from measurements in similar examinations using other imaging methods in the same or different subjects.

In contrast to the limited spatial information, EEG and MEG have a very high resolution in time. This unique feature of the EEG and MEG to measure rapid transients and oscillations of brain activity allows for the imaging of source current signals or source waveforms over time (4D images). The repeated measurements at different time instances (epochs) also provide additional information which can be used to decrease the voxel size or to localize a larger number of discrete active voxels.

The iterative procedure may be made even faster and more stable by starting with a coarse mesh and decreasing the voxel size during the procedure. Hence, one may start by dividing the brain into a mesh of large voxels, identifying the voxels containing active regions, and then iteratively subdividing the remaining active voxels into smaller voxels before repeating the iterative core loop of the discrete source scan procedure (Fig. 2).

After this procedure has terminated, the second method of the invention, the 'Multiple Source Image Procedure' (Fig. 3), can be used to create a voxel image with the magnitude of the source activity at each time instance or with the magnitude and phase of the source activity at different frequencies. The spatial 3D images are determined by scanning the whole mesh using a so-called "probe source" at each voxel in conjunction with the sources at the centers of the 'active' voxels. Using the probe source scan on a fine image mesh, the position, orientation and strength of the source activity in the voxels identified as active voxels may be optimized. The estimation of the source activities can be improved by analyzing the contents of the continuously measured detector data in the time and frequency domains. In addition, if some knowledge exists on the temporal-spectral features of the source activities, spectral-temporal filters can be constructed that increase the source signals while reducing or suppressing the noise signals that are also picked up by the detectors. These noise signals may originate inside or outside of the 3D volume or can arise in the physical detector and/or the measurement apparatus.

Therefore, sharp discrete events in time can be specified that are associated with the emission of a source signal from a discrete location within the 3D volume. Such emitted source signals are generally limited in duration and have been named, for example, 'atoms' of the EEG or MEG (Miwakeichi et. al. 2004).

Any emitted source signal can be decomposed into a weighted sum of basic patterns defined by a) a 2D set of complex wavelet functions specified by the 2 dimensions frequency and duration, or the number of cycles at that frequency, respectively. Each wavelet function is a complex signal function having a cosine and a sine part with an envelope of a specific duration. b) a 2D set of predefined patterns specified by duration and the numbers of underlying periodic oscillations. These patterns are derived either from knowledge of previous measurements and analysis from other subjects, or they are determined using the reverse spectral transform of the spectral eigenvectors of covariance matrices of the source signals calculated progressively over longer epochs A, B etc. of a long-term detector measurement (cf. Figs. 4 & 8).

Therefore, the third method of the invention, Time & Pattern Definition Procedure' (Fig. 4), can be used to estimate the discrete time of occurence, the duration and the spectral content or periodicity parameters defined as the peaks of the weighting function of the basic patterns for each source activity. The weighting function determines the amplitudes of each basic pattern in the above defined sets of predetermined patterns in a linear combination of the basic pattern signals that constitute the source signal. The parameters time of occurence, duration and spectral content or periodicity, respectively, as well as the weighting function can be determined using the optimization procedure defined below.

Finally, source activity events can be defined by a small set of parameters using both the spatial parameters location and orientatior of each source defined by the discrete source scan in conjunction with the temporal, duration and frequency or periodicity parameters, respectively, obtained by the time & pattern definition procedure in a combined maximum likelihood optimization procedure defined below. Thus, the detector data can be optimally filtered in space, time and frequency to estimate more robust source locations, orientations and source signals for each epoch and sub-epoch while largely reducing the influence of noise.

Detailed Description of the Drawings and the Mathematical Background

Fig. 5 illustrates the physical relationship between n detectors on a 2-dimensional surface and m sources within a 3D physical volume. The detectors at locations x_k (k=1...n) over or in the 3D volume pick up time-varying data dκ(t) that result from the propagation of source activities s,(t) of m different, circumscribed or discrete source regions (compartments or voxels) with center locations r, (i=1...m) within the 3D volume to the n detectors. The number m, locations r,, and the time-varying activities s,(t) of the source regions are unknown. A solution of this discrete linear inverse problem is obtained under the following assumptions: a) A model of the physical transfer function T(k,i) is available that relates the strength of the source signal s,(t) at each voxel i (i=1...m) to the strength of the transmitted signal at each detector k (k=1...n). b) The number of discrete sources is less than the number of detectors times the number of statistically independent observations over time. A conservative, robust estimation is achieved if the number of simultaneously active sources to be estimated is kept less than the number of detectors. c) The signals dk(t) at the detectors are the linear sum of the attenuated source signals s,(t) according to the laws of physics and their sum is non-zero at each detector (Fig. 5):

(1 ) d_k(t) = ∑,=i. _m T(k.i) s,(t) + n_k(t)

In general, each detector will also measure some noise signal described by nι_<(t).

If the observations are made at discrete times t = t|....t_Nt, then this superposition equation can be formulated as a matrix equation (in the following, bold lower case letters designate vectors, e.g. d(t), and bold upper case letters designate matrices, e.g. D. All non-bold letters are scalars): (2a) D = T S + N (i.e. D_M = ∑,=i, m T _k|l S ,,_t + N_k,_t) In vector annotation this can be written as: (2b) d(t) = T s(t) + n(t)

D is a matrix with n rows and N_t columns (Fig. 5). The rows represent the time-varying signals of each detector (k=1... n) and the columns are the observation vectors at one time sample with N_t discrete time instances (t = t, ....t_Nt).

S is a matrix with m rows and N_t columns. The rows represent the unknown time-varying signals of each source (i=1...m) and the N_t columns are the source activity vectors at one time representing N_t discrete time instances (t = ti....t_Nt).

T is the transfer matrix or forward model matrix with n rows and m columns. Column i represents the topography of source i at the surface, i.e. the transfer coefficients at each detector for a source activity of unit strength at location r,. Therefore, T contains the attenuation factors for the propagation of unit activity from each source i to each detector k (Fig 5.)-

For any T calculated for m sources in different compartments or voxels (Fig. 5), T has rank m and the source activities s(t) can be calculated by operating onto equation (2a) from the left with the linear inverse operator T^"1 defined by

(3a) T^"1 = (T T) ^"1 T with T¹ designating the transpose of T.

To reduce the effect of noise in the probe source images (Fig. 3), a smoothing or regularization vector p can be inserted prior to calculating the inverse of the inner product of T:

(3b) T¹ = (T¹ T + p l) ^"1 T"

The source activities can be estimated using an implicit least-squares minimization of the noise term by

(4a) S = T^"1 D in matrix annotation, or

(4b) s(t) = T^"1 d(t) in vector annotation

Source activity may arise in the 3D volume spontaneously (Figs. 5 & 7) or following an input signal (Fig. 6) that represents, for example, an external stimulus, or time-locked to an output signal, e.g. a response (Fig. 6). Accordingly, the detector measurements over time can be divided at least into two epochs, e.g. epoch A and epoch B, such that epoch B represents some baseline epoch prior to the spontaneous or event-related activity and A the epoch when these activities occur within the 3D volume. Epoch C can be defined as the response-locked epoch, for example.

To find the unknown source locations, let the volume be divided into a finite number of discrete compartments. The mesh may consist of regular cubic voxels (Fig. 5), tetrahedrons, as used, for example, for finite element models (FEM) of the volume conduction in the head or heart, or more complex polyhedrons.

First, consider the situation that 3 discrete sources contribute to the detector measurements and that each equivalent source describes the local source activity in the surrounding volume completely. In this case, an additional probe source p (Fig. 5, black circle) that is placed at any other voxel will not record any source activity of sources 1-3 since the inner product of the inverse and forward vectors of T yields 1 for the same source (i = j) and 0 for different sources (i ≠j)

j

In other words, if sources are found at the correct equivalent locations, i.e. at the centers of the source regions, then there will be no cross-talk from one source to another, and a randomly placed probe source will only pick up noise and no source signal from the other sources.

/ Measurement of the Estimated Source Activity Value (Q-value)

To obtain a measure of estimated source activity, hereafter called Q-value, for each source voxel i, one can either use the maximum over time of the estimated power of the source activity during epochs A and B

(6a) Q_Ap² = maximum(t) (SAI² (t) , t=1 N₁ ) = s_A,² (t_A-max), and

(6b) Q_Bp² = maximum(t) (s_Bl ² (t) , t=1 N₄ ) = s_B,² (t_B-max), or one can use the mean power P over time (6c) Q_Amp² = mean(t) (s_Al ² (t) , t=1 ,....Nt ) = < s_Al ² >, and (6d) Q_Bmp² = mean(t) (s_Bl ² (t) , t=1 ,....Nt ) = < s_Bl ² > where <s_A,²> denotes the mean of s_A ² over epoch A, etc.

The mean power over time can be calculated effectively for each data set A and B by using the spatial covariance matrix C(t) over time in the detector space with dimension n. C(t) is given by the inner product of the data matrix D with its transpose (7) C = D D' = U ΣV V Σ^I U^I = U Σ Σ U^I

By applying a singular value decomposition to D, as in eq. 7, and splitting the signals of the n-dimensional detector space into a signal subspace with N_s eigenvalues (first N_s entries of the diagonal matrix ∑_s above a threshold tr_ev, and into a noise subspace with N_N eigenvalues (n = Ns + N_N) according to Mosher et al. (1992), D can be separated into

(8) D = U_S∑S V_S + UN∑NVN

The mean source power of the signal and of the noise can then be computed by projecting the inverse of T onto the signal and noise subspaces:

(9a) <ss,²> = trace, (T¹Us ∑s ∑s' U_s' T¹') (9b) <S_N,²> = trace, (T¹U_N Σ_N Σ_N' U_N' T¹')

When comparing two epochs A and B, for example baseline B and signal epoch A, or the signal and noise part of the data matrix D defined by the eigenvectors Us and U_N, it is convenient to calculate a ratio of the estimated activity between both data sets (setting S=A and N=B):

(1 Oa) Q² _pr (i) = S_A,²(t A-max) / <S_B ,²>

(10b) Q²pr (i) = <SA,²> / <SB ²>

As an alternative to a baseline epoch B prior to a signal epoch A, either the signal subspace spanned by U_N ∑N or the plus-minus average (Schimmel 1967) can be used to define the signals B, if a sufficient number of event-related sub-epochs A is available before averaging or if sub-averages have been computed that can be subtracted to estimate the noise signals B.

Then, one can define an amplitude ratio of the signal to noise or an amplitude ratio to compare two epochs A and B that might reflect two different task conditions during the measurement, similar to blocked designs in fMRI taking the square root (sqrt) of eq. 10:

(11a) Q(r,) = sqrt (<S²AI>) / sqrt (<s² _Bl>) -1 , if the mean activities over epochs or conditions A and B are considered and if <SA²> > <s_B,²> , or

(11 b) Q(r,) = sqrt (<s² _Bl>) / sqrt (<s² _Al>) -1 , if the mean activities over epochs / conditions A and B are considered and if <SB ²> > <s_Aι²> , or

(11c) Q(r.) = |s_A,|(tA-maχ) / sqrt (<s² _Bl>) -1 , if the peak in SA is considered and if s_Aι²(tA-max) > <SB_I ^2>, i.e. the peak source activity in epoch A is larger than the mean source activity in epoch B, or,

(11 d) Q(r,) = -|s_B,|(t_B-maχ) / sqrt (<s² _Al>) -1 , if the peak in SB is considered and if Se^ts-max) > ^<SAI^2>, i.e. the peak source activity in epoch B is larger than the mean source activity in epoch A.

These Q-values are equivalent to a z-score -1 , i.e. the source signal is divided by a standard deviation and 1 is subtracted.

It can be convenient to define the Q-value on the basis of this amplitude ratio using a logarithmic scale in decibels (= 20 log (x)) by taking the root (= implicit division by 2 outside of the log (x)) and the logarithm of the source power. Thus, one can rewrite eq. 10 as (again setting S=A and N=B):

(11e) Q(r.) = 10 log <s² _A,> - 10 log <s² _Bl> , if the mean activities over epochs or conditions A and B are considered, or, (11 f) Q(r,) = 10 log s_A,² (t_A-max) - 10 log <s_B,²> , if the peak in SA is considered and if s_A ²(tA-max) > <SBI^2>, i.e. the peak source activity in epoch A is larger than the mean source activity in epoch B, or,

(11g) Q(r,) = 10 log <s_Al ²> - 10 log s_B,² (t _B-maχ). if the peak in SB is considered and if SB²(tB-max) > ^<SA^2>, i.e. the peak source activity in epoch B is larger than the mean source activity in epoch A.

For the search of local maxima in Q(r,) one can use the Laplacian operator in 3D: (12a) Lap-3D (Q) = d² Q / dx² + d² Q / dy² + d² Q / dz² ,

To estimate local maxima of the source activity s,(t) in space and time, one can use the Laplacian operator in 4D:

(12b) Lap-4D (s,(t)) = d² s,(t) / dx² + d² s,(t) / dy² + d² s,(t) / dz² + d² s,(t) / dt²

To obtain estimates of the source signals for the final spatio-temporal 4D images, i.e. time-varying 3D images that are less contaminated by the noise, the source activity matrix S may be calculated using an inverse of the transfer matrix T that is projected into the 3D signal space spanned by the signal eigenvectors or topographies U_s of D (Scherg and Berg, 1996):

(13) S = (T¹ T)-¹ T¹ Us Us¹ D Discrete Source Scan Procedure

As a consequence of the inverse principle specified by eq. 5 and in Scherg (1992), the first method of the invention defines the 'Discrete Source Scan Procedure' (Fig. 2) to search the finite voxel set for active and inactive voxels such that the active voxels contribute signal to the detectors whereas the inactive voxels contribute mainly noise, i.e. signals that are considerably smaller in magnitude as compared to the signals of the active voxels. The discrete search process over a finite mesh yields a unique result, if the above assumption b) is valid.

Principally, there are three different methods to determine the active voxels. These methods may also be used in combination by

A. using a function that minimizes the difference between the measured detector data and the predicted or modeled data contributed by the active voxels, i.e.

(14) rv =|| D - T S ||^p = || D - T T^"1 D ||^p = || ( I - T T¹ ) D ||^p where rv is the residual, unexplained variance, I is the unit matrix and p is the order of the norm that is typically set to 2 (least-squares minimization) or 1 (absolute magnitude minimization). Since the transfer matrix T = T (rj,0j,øj) depends non-linearly on the locations η of the m active sources (and on the orientation angles θ, and φ_\ in the case of single dipole vectors) these spatial parameters are traditionally determined by an iterative procedure that minimizes (14) using downhill or simplex fitting algorithms, genetic algorithms or simulated annealing (Scherg 1990; Mosher et al. 1992; Scherg and Berg 1996).

In the discrete source scan method of this invention, a related, but new and different optimization function is used for the local fine-tuning of the source parameters r_lt θ_\, and φ_\ after active voxels have been identified. This optimization is based on the following maximum likelihood (ML) function consisting of 5 terms or fit criteria:

(15) ML = w_Q |Q| - W_n, In (rv / ||D D'||^p + COnSt_n,) - w_s ln ( ||S S'||^P / ||D D'||^p + const_s)

⁺ In w_md ( minimum ( min-distancβi, [ ( Tj ₁ Ti ) ; const_md ) )

+ In w_sc ( minimum ( Σ (all i ≠\ ) (1- tⁿi't^πι ) ; const_sc ) )

The weights W_x of each of the fit criteria can be selected to emphasize one or more criteria in the fitting procedure. The constants const_x serve to reduce the effects of criteria that have already been optimized. The Q-value is always optimized. The source parameters r, θ_u <p_\ are locally optimized using the parameters of the m active voxels determined by method B or C (below) and the ML function (15). The optimization uses an iterative procedure based on a) the simplex fitting algorithm with small local starting step size, or b) on a genetic algorithm adapted for local optimization, or c) on simulated annealing.

B. iteratively using a combination of q randomly selected voxels and p=m-q selected sources (Fig. 2) that were defined a) as local maxima during previous iterations, or b) by systematically scanning the voxels of the mesh, or c) by restraining the random selection to those voxels that belong to the last list of active voxels. During this iteration procedure a Q-value is calculated for each voxel in the set as defined above by eq. 11 , or alternatively, as defined below in C. Furthermore, storages of the maximum, minimum and mean Q-values are kept and updated at specified intervals or under specific conditions, for example, if rv or ML exceed predefined or adaptive thresholds tr™ or tr_Mι_- The storages are evaluated at these intervals or under these conditions to update the lists of active and inactive voxels. The iterative discrete source scan procedure terminates if a) the number of active voxels m is less than the number of detectors n, b) a certain number of iterations has elapsed, or c) a number p of selected voxels with peaks, i.e. local maxima that are separated by a minimum spatial distance threshold t_md is less than, for example 2S, i.e. twice the number of dimension of the signal space defined by Us- The threshold 2S is reasonable in the case of EEG and MEG measurements because of the high likelyhood of symmetric and nearly-synchronous activities in the right and left hemispheres.

C. using a 'Multiple Source Beam Former' (MSBF) in an iterative procedure that consists of the following steps: a) The starting iteration 1 uses a dual beamformer in order to detect local maxima of the Q-value across a predefined mesh of compartments or voxels (Fig. 2). For the analysis of EEG and/or MEG data it is useful to start iteration 1 with a bilaterally symmetric beamformer (BSBF) to match the similarities of the right and left hemispheres of the brain. b) A regional source is placed into each local maximum of the Q-values (or its Laplacian according to eq. 12) that has been determined during the previous iteration and that is above a predefined or adaptive threshold tr_Q , tr_Lap-3D(Q) tr_{Lap- 4D(S}i₍D₎ according to eqs. 11-12. c) Optionally, the regional source is converted into a single dipole by rotating the 3 vectors of the regional source such that the first vector explains the maximum of the projection into the signal space Us or such that it explains the source activity s_Al ²(t A/B-maχ) completely at t = t A/B-max - d) During each subsequent iteration, the sources determined at the local maxima in steps b) and c) are kept as active sources in conjunction with an additional probe source that is used to perform the MSBF scan. The scan consists of the calculation of a Q-value for the probe source(s) at each compartment or voxel according to eq. 19 below. e) For local fine-tuning of the maxima, the location of each maximum stored for the next iteration may be locally randomized within a radius η_oc around the local maximum. f) The predefined thresholds may be adapted. g) The iterative procedure runs repeatedly over steps b), c), d), e) and f) and stops if

1 ) no more maxima are detected above the predefined or adaptive thresholds during a predetermined number of successive iteration, or

2) the number of determined sources m is equal to the number of detectors n, or

3) the number of determined sources m is equal to or greater than 2S, i.e. twice the number of dimensions of the signal space defined by Us, or

4) a predetermined number of iterations has been reached.

The inverse spatial operator W¹ of the single source beamformer yields the following spatial filter output when operating on the forward topography transfer matrix T (van Veen 1997, eq. 12):

(16) W¹(r_q) T(r,) = I , if r_q = r, , & W^"1(r_q) T(r) = 0 , if r_q != r,

As a consequence of this condition, the correlations over time or frequency between the activities of the sources r_q and r, as defined by the spatial covariance matrices C(t) and C(f) are required to be zero. If the sources are correlated, then the single source beamformer identifies maxima at erroneous locations. The amount of the localization error depends on the degree of correlation. The method C of the discrete source scan procedure of this invention avoids the requirement of non-correlatedness by using multiple sources in a beam former operator modified to accommodate multiple sources. While van Veen's inverse operator W(r_q) (US 5,263,488; van Veen 1997, eq. 23) is a spatial filter corresponding in function to the inverse operator T^"1 in the previously outlined methods, it is not identical to the T^"1 defined by one regional source in eqs. 3 & 13:

(17) W(r_q) = [ T(r_q) σ¹ T(r_q) ] -¹ r(r_q) σ¹

Furthermore, the T transfer matrix of van Veen consists only of the 3 forward topography vectors of one regional source. In contrast, the method C of the discrete source scan procedure of this invention extends van Veen's formulation by specifying a transfer matrix T that uses the topographies of multiple dipoles and/or regional sources in conjunction to form a multiple source transfer matrix T and to compute a multiple source beamformer (MSBF) by

(18a) Q_p ²(n) = trace, { [ T¹ C-¹ T ] -¹ J summing only those elements in the trace which belong to source i, or

(18b) <ss,²> = trace, [ T U_s ∑s^"1 ∑s^"1' U_s' T ] ^"1 calculating the mean source power of the signal space, similar to eq. 9, or

(18b) <s_Nl ²> = trace, [ T¹ U_N Σ_N ^'1 ∑N^"1' U_N ¹ T ] ^"1 calculating the mean source power of the noise space, similar to eq. 9, ,or

(18c) s_A,²(t_A-max) = H [ T¹ UA ∑A^"1 ∑A^"1' UA¹ T¹ ] ^"1 T¹ U_A ∑_A ^"1 d_A(t_A-max) Il calculating the source power in epoch A at its peak time t_A-max-

Similarly to eqs. 10 and 11 let S=A and N=B to calculate the Q-value using the MSBF on the basis of the source activities defined in eq. 18:

(19a) Q(r,) = sqrt (<S²A,>) / sqrt (<s² _Bl>) -1 , if the mean activities over epochs / conditions A and B are considered and if <SA ²> > <s_B ,²> , or

(19b) Q(r,) = sqrt (<s² _Bl>) / sqrt (<s² _Al>) -1 , if the mean activities over epochs / conditions A and B are considered and if <s_Bl ²> > <SA,²> , or

(19c) Q(r,) = |s_Aι|(t_A-max) / sqrt (<s² _Bl>) -1 , if the peak in s_A is considered and if s_Al ²(tA-_maχ) > <SB ^2>, i.e. the peak source activity in epoch A is larger than the mean source activity in epoch B, or,

(19d) Q(r,) = -|s_Bl|(t_B-max) / sqrt (<s² _Al>) -1 , if the peak in s_B is considered and if s_Bι²(tB-max) > <SA_I ^2>, i.e. the peak source activity in epoch B is larger than the mean source activity in epoch A.

Alternatively and according to eq. 11 , the Q-value can be defined on the basis of this amplitude ratio using a logarithmic scale in decibels:

(19e) Q(r,) = 10 log <s² _Al> - 10 log <s² _Bl> , if the mean activities over epochs or conditions A and B are considered, or, (19f) Q(r.) = 10 log s_Al ² (t_A-max) - 10 log <s_Bl ²> , if the peak in s_A is considered and if s_A ²(tA_-maχ) > <s_B ²>, i.e. the peak source activity in epoch A is larger than the mean source activity in epoch B, or,

(19g) Q(r,) = 10 log <s_Al ²> - 10 log s_Bl ² (t _B-max), if the peak in s_B is considered and if s_Bι²(t_B-rnaχ) > <s_Aι²>, i.e. the peak source activity in epoch B is larger than the mean source activity in epoch A.

Furthermore, the Lap-3D and Lap4D operators can be defined for the Q-value of the MSBF according to eq. 12.

In summary, the method C of the discrete source scan procedure of this invention extends van Veen's beam former from one regional source to several dipoles and/or regional sources, thus dropping the intrinsic requirement of non- correlatedness of the source activities over time. Similarly, it extends the SAM method (WO00/10454) that uses only one oriented dipole source at each grid point to identify uncorrelated sources of the MEG to a method utilizing a higher number of sources that can be correlated.

Multiple Source Image Procedure

The second method of the invention, the 'Multiple Source Image Procedure' can be used to create a voxel image with the magnitude of the source activity at each time instance or with the magnitude and phase of the source activity at different frequencies. As shown in Fig. 3, the inputs to the image procedure are: a) spatial data vectors either directly from the data matrix D or from spatial eigenvectors U of covariance matrices, b) a source model with regional sources at locations η or dipoles at locations r with orientations o, where Oj is the unit orientation vector in Cartesian coordinates.

The spatial 3D images are determined by scanning a fine mesh that covers the whole 3D volume using a so-called 'probe source', i.e. either a regional source or a dipole with fixed orientation, in conjunction with the other sources at locations η. The Q-value or Laplacian of the Q-value of the probe source is evaluated according to eqs. 11 or 12 and eq. 3b at each compartment or voxel of the mesh, scaled and displayed as a 3D image.

When using the Q-value output of the MSBF according to eq. 19, a combined 3D image is created from m 3D images that are calculated for each of the m sources by excluding the respective source and adding the probe source such that the maximum of the Q-value over the m images is selected for each voxel in the combined 3D image.

Source and Detector Data in the Frequency and Time-Frequency Domains

Both the detector and the source signals of one or more time epochs can be transformed into the frequency domain. This is achieved by multiplying eqs. 2 & 4 from the right with a time to frequency transformation operator F, defined e.g. by the coefficients of a fast Fourier (FFT) or a wavelet transform:

(20a) S F = T¹ D F with F being a matrix with N_t rows and N_F columns. NF <N_t /2 is the number of basic frequencies or wavelets used in the transform. Each element of F is a complex number described by a cosine (real) and sine (imaginary) part. If desired, the number of time samples N_t can be reduced prior to applying F by low-pass filtering, and - in the case of event-related inputs - taking the logarithm of time, and decimating the temporal signals contained in D. Due to the linearity of the transformation eq. 20a can be rewritten as:

(20b) SF = T^"1 D_F in matix annotation, where D_F and SF represent the complex data and source activity matrices, respectively, in the frequency domain, or:

(20c) S_F(f, Φ) = T^"1 dF (f, Φ) in vector annotation, where S_F(f, Φ) is defined as a function of frequency f and phase φ.

If a moving spectral representation is created over an epoch with several sub-epochs, D_F and SF become functions of time and describe the temporal spectral evolution (TSE) over time: (2Od) S_F(t) = T¹ DF (t) in matrix annotation, or

(20c) s_F(t, f, ψ) = T^"1 d_F (t, f, φ) in vector annotation, i.e. the data and source activities are represented as functions of time t, frequency f, and phase φ. The spatial covahance matrices in the frequency and time-frequency domains then become:

(21 ) C(f) = D_F D_F*

(22) C(f, t) = D_F (t) D_F* (t)

The elements of the matrices D_F , S_F , C(f) and C(f, t) are complex values with real and imaginary parts and * denotes the conjugate complex. The main diagonal of C_F is real valued since C_F is a Hermitian matrix.

Using the same formulations as outlined above, the power, magnitude, phase, and other derived parameters of the source activity matrix S_F(t) can be computed and displayed using time-frequency plots. For example, the log (magnitude) of source activity is defined in decibels as:

(23) |s,(f,t)| [dB] = 10 log [ (V¹ d_F(f, t) ) , (tT¹ d_F(f, t))^*, ]

Time & Pattern Definition Procedure

The estimation of the source activities can be improved by analyzing the contents of the continuously measured detector data in the time and frequency domains. As outlined in Figs. 4 & 8, temporal and spectral vectors, or eigenvectors of the various covariance matrices can be prepared that represent the underlying temporal and spectral patterns of the measured data. The data are sub-divided into overlapping segments of short duration, for example 2 second segments that are evaluated in steps of 1 second. The mean of each segment is subtracted and the resulting signals are tapered using a cos² function centered on the segment with a period corresponding to the segment length. The tapered signal of each detector k or estimated source i is then transformed into the frequency domain using a Fourier or a wavelet transformation. This transformation has to be applied only once to the detector data, since the source activity matrix S_F in the frequency domain can be computed according to eg. 20.

Similar to the definition of the spatial covariance matrices in eqs. 7 & 21 , temporal covariance matrices C_t of dimension N_t • N_t can be defined by the outer product of the temporal detector or source signal vectors of each sub-epoch A_p (or B_q) having N_t elements:

(24a) Ct (d_kAp) = dkAp¹ d_kAp

(24b) Ct (s,Ap) = S₁Ap¹ s_kAp = D_Ap" V¹' t^'1D_Ap

Similarly, spectral covahance matrices of dimension NF ^* N_F , N_F <N_t / 2 can be defined for each sub-epoch A_p after tapering and Fourier or wavelet transformation:

(25a) CF (dAkp) = F¹ c-Akp' dAkp F = CJFAIΦ¹ dFAkp

(25b) CF (SAIP) = F¹ s_Aιp' s_Aιp F = SFAIP' S_FAI_P = D_FA_P' t,^"1' tr¹D_FAp

The complex spectral covariance matrices are averaged over the N_p sub-epochs of condition A and the average is decomposed into complex spectral eigenvectors, or basic spectral patterns VFIA-

(26) CF (SAI) = 1 / Np ∑p CF (SAIP) = VFAIP ∑FAIP ∑FAIP VFAIP*

Eq. (26) can similarly be computed for each detector channel k or averaged over all detector channels.

Finally, the reverse Fourier transform of each spectral eigenvector V_FIA_J of dimension N_F defines a basic pattern in time that is specific for the signals emitted from source i:

(27) v_tlJ = F^"1 v_FAlJ

A specific subset or linear combination v_F,h of the N_F eigenvectors v_Fιj of any epoch A, B or C can be used to create an inverse spectral filter to transform the spectral vectors d_F and/or s_F dividing by the specific pattern v_FAιh . and taking the reverse Fourier or wavelet transform into optimally filtered time signals:

(28a) d_t = F^"1 ( d_F / v_FA,h ) and,

)

Using the decomposition of the overlapping segments or sub-epochs and their transformation into the frequency domain as described above, the inverse filter of eq. 28 can be applied to each detector and source signal in order to recompose optimally filtered detector data over time. The resulting filtered signals will then exhibit sharp transients, similar to delta functions, whenever an event occurs that is described by a specific linear combination of the basic spectral patterns v_FAιh- These transients can then be used together with appropriate threshold criteria to detect, for example, spikes or seizure onsets, or slow waves in abnormal EEG and/or MEG data recorded from patients with epilespy, stroke, migraine and/or other disorders.

The basic spectral pattern matrix V_flA of source i forms an orthonormal base to describe the measured source signals during condition A. Alternative orthonormal pattern bases can be defined directly using the basic sine and cosine functions of the Fourier or wavelet transforms, or a specific set of patterns V_t that form a basis for the waveforms of all sources. Assume that there are N_b patterns and their spectral transform is

(28) V_F = V_t F with V_F having Nb rows and N_F columns and V_t having N_b rows and N_t columns. The source signals in the frequency domain can then be described as a product of a discrete time-shift operator i transformed to the frequency domain p_F, = p, F with a weighting function that specifies the combination of basic patterns for each source, i.e.

(29a) S_F1 = P_R V_F ¹ W₁ ¹ and after reverse Fourier transform into the time domain (29b) s, = F^"1 p_Fl V_F ' w,¹

The time-shift operator p, consists of a time vector delta function δ (i-τ) of length N_t with zero entries except for a unit entry at the time t=r, defined by the onset or by the first peak of source activity i:

(30) p,(t) = δ(t-_Tι)

Thus, the source activity time vector s, can be constructed using a discrete time of occurrence η and a weighting vector of basic patterns w,.

The basic set of patterns does not necessarily have to consist of orthonormal basic function. It can be considered, for example, as a 2D array of basic waveforms that are arranged and sorted according to duration and perodiocity or signal frequency, respectively, such that the maxima in the weighting vectors can be interpreted as the duration and periodicity or frequency parameters of each source activity event that reflects a time-limited emitted pattern.

The times of occurrence r, and the weighting vectors w, can be estimated in an interactive least squares (or ML) fitting procedure in analogy to eqs. 14-15 minimizing:

(31 ) rv =|| D - T S ||^p = || D - D S^"1 S ||^p = || D ( I - S^"1 S) ||^p where S(7i,w,) is a function of the times of occurrence r, and the weighting vectors w, and taking into account that an estimated transfer matrix T_e can be calculated by inverting eq. 2 from the right:

(32) T_e = D S^"1

Finally, a mutual optimization procedure can be initiated that optimizes both the estimated transfer matrix T_e (as computed by eq. 32) and the source activity matrix S_e (as computed by eq. 2) by varying or fitting the parameters of the activity events that define the forward model transfer matrix T_m and the model-constructed source activity S_m. This is achieved by optimizing the following maximum likelihood criterion and using equations 14 and 31 :

(33) ML = - In H D ( I - S_m ¹ S_m) ||^p - In || ( I - T_m T_m ^"1 ) D ||^p

+ In Σ, te, W (IIt₈₁H || t_mi||) + In Σ, s_eι s_mι / (||s_e,|| || s_m,||)

Source Coherence

For more generality, the time-frequency annotation g(f,t) and C(f,t) is used in the following for all signals g and matrices C.

The coherence coh²(r_ι,r,,f,t) between 2 sources i and j can be obtained by projecting the spatial complex covariance matrices C(f,t) for each frequency f and time epoch t onto the inverse spatial operators V¹ and V¹ (cf. Gross et al. 2001 ):

(34) c_Fs,(f ,t) = tr¹ C(f,t) t, -^{1 t}, where c is a complex number described, e.g. by amplitude and phase c= |c| e *^ω+φ , and calculating

(35) coh^.r_j.f ,t) = (c_sι(f,t) c_s,(f ,t)^*)² / c_sl(f ,t)² c_SJ(f ,t)²

Coherence matrices can be summed in the time-frequency domain over 'a' discrete frequencies and 'b' times to obtain more information on multiple sources for the discrete source scan procedure

(36) coh^.r_j.frUrt_b) = ∑_pq (c_sl(f_p,t_q) c_s,(fp,t_q)^* f I c_s,(f_p,t)² c_SJ(f_p,t)²

The coherence of a source signal s,(f,t), with an external signal g(f,t), another detector signal g(f,t)= dι(f,t), or another source signal g(f,t)= s/f.t), is defined by the correlation in the frequency domain. For each sub-epoch p of the measured data the coherence term is defined as: (37) cohp = (d_kp(f,t) g(f ,t)* f I [ (d_kp(f ,t) d_kp(f,t)^* ) ( g_p(f ,t) g_p(f,t)* ) ] where g(f,t)* is the conjugate complex of g(f,t).

Normally, coherence is summated over p sub-epochs to yield the correlation of epochs in the frequency domain. Thus, the source coherence of source signal S_Fi(f.t) with g_F(f,t) in the time-frequency domain is calculated by:

(38) coh(s_R(f,t), g_F(f,t)) = Σ_p(s_Fip(f,t) g_Fp(f,t)^* f I ∑_p s_Fi(f,t) s_Fi(f,t)^* ∑_P g_Fp(f.t) g_Fp(f,t)^* =

= ∑p ( ti-¹D_Fp(f, t) g_Fp (f,t)^* f I ∑_P ti-¹D_Fp(f, t) (ti-¹D_Fp(f, t))^* ∑_P g_Fp(f,t) g_Fp(f,t)^*

An examplary hardware implementation of the present invention is schematically illustrated in Fig. 11. A measurement unit 20 picks up the measurement data from the plurality of detectors 10, e.g. EEG detectors arranged around the scalp of a patient. The data are then fed to a processing unit 30 to convert the analog signals to digital numbers, to store the converted digital data and to perform the calculation and analysis steps. The processing unit comprises an analog-to-digital ASC converter unit 31 , a storage unit 32 and a calculation unit 33. The results are then forwarded by an interface unit 34 to a display unit 35 for display and/or interactive analysis of the results (cf. Fig. 1 ).

References

Dale A.M., Liu, A.K., Fischl, B. R., Buckner, R.L. Belliveau, J.W., Lewine, J. D., Halgren E. (2000). Dynamic statistical parametric mapping: combining fMRI and MEG for high-resolution imaging of cortical activity. Neuron, 26: 55-67.

Darvas, F., Schmitt, U., Louis, A.K., Fuchs, M., Knoll, G., Buchner, H. (2001 ). Spatio-temporal current density reconstruction (stCDR) from EEG/MEG data. Brain Topography 13: 195-207.

Gross, J., Kujala, J., Hamalainen, M., Timmermann, L., Schnitzler, A., Salmelin, R. (2001 ). Dynamic imaging of coherent sources: studying neural interactions in the human brain. PNAS 98: 694-699.

Hoechstetter, K., Bornfleth, H., Weckesser, D., IHe, N., Berg, P., Scherg, M. (2004). BESA source coherence: a new method to study cortical oscillatory coupling. Brain Topography 16: 233- 238.

Liu, A., Dale, A.M., Belliveau, J.W. (2002) Monte Carlo simulation studies of EEG and MEG localization accuracy. Human Brain Mapping 16: 47-62. Miwakeichi, F., Martinez-Montes, E., Valdes-Sosa, PA, Nishiyama, N., Mizuhara, H., Yamaguchi, Y. (2004). Decomposing EEG data into space-time-frequency components using parallel factor analysis. Neurolmage 22: 1035-1045.

Mosher, J. C, Lewis, P. S., Leahy, R.M.(1992). Multiple dipole modeling and localization from spatio-temporal MEG data. IEEE Trans Biomed 39: 541-557.

Scherg, M. (1984). Spatio-temporal modelling of early auditory evoked potentials. Rev. Laryngol. (Bordeaux), 105: 163-170.

Scherg, M. and von Cramon, D. (1985). Two bilateral sources of the late AEP as identified by a spatio-temporal dipole model. Electroenceph. din. Neurophysiol., 62: 32-44.

Scherg, M. and von Cramon, D. (1986). Evoked dipole source potentials of the human auditory cortex. Electroenceph. Clin. Neurophysiol. 65: 344-360.

Scherg, M. (1990). Fundamentals of dipole source potential analysis. In: Auditory evoked magnetic fields and electric potentials (eds. F. Grandori, M. Hoke and G. L. Romani). Advances in Audiology, Vol. 6. Karger, Basel, pp 40-69.

Scherg, M. and Picton, T. W. (1991). Separation and identification of event-related . potential components by brain electric source analysis. In: Event-related potentials of the brain (eds. C. H. M. Brunia, G. Mulder and M.N. Verbaten). EEG Suppl. 42, Elsevier Science B.V., Amsterdam, pp 24-37.

Scherg, M. (1992). Functional imaging and localization of electromagnetic brain activity. Brain Topography 5: 103-111.

Scherg, M. and J. S. Ebersole (1993). Models of brain sources. Brain Topography 5: 419

Scherg, M. and Berg, P. (1996). New concepts of brain source imaging and localization. In: Functional Neuroscience (eds. C. Barber, G. Celesia, G. C. Comi and F. Mauguiere). Elsevier Science B. V., Amsterdam, EEG Suppl. 46, pp 127-137.

Scherg, M., Bast, T., Hoechstetter, K., IHe, N., Weckesser, D., Bornfleth, H., Berg, P. (2004). Brain source montages improve the non-invasive diagnosis in epilepsy. International Congress Series, 1270C: 15-19.

Van Veen, B.D., van Drongelen, W., Yuchtman, M., Suzuki, A., (1997). Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans Biomed 44: 867-880.

Claims

Claims:

1. A method of determining a finite number of discrete sources within a 3- dimensional volume, comprising the steps of: a) measuring a physical quantity using a number of n detectors having a known position, b) dividing the 3-dimensional volume into a plurality of discrete compartments, c) selecting a subset of m compartments, d) determining estimated activity values (Q-values) of the m compartments based on the measurement data, e) storing and updating the estimated activity values, f) iteratively repeating steps c), d) and e) with new subsets of compartments, g) determining a number of compartments as active sources based on the stored and updated estimated activity values.

2. The method of claim 1 further comprising the step of calculating the source activity of the determined active sources using the measurement results of the n detectors.

3. The method of claim 2, further comprising the step of generating an image of the scanned volume based on the activity values of the active sources.

4. The method of one of claims 1 to 3, wherein the subset of m compartments is selected randomly in step c).

5. The method of one of claims 1 to 4, comprising utilizing prior knowledge about expected activity within the volume for selection of the subset of m compartments in step c).

6. The method of one of claims 1 to 5, wherein step g) includes determining a compartment as an active source if the estimated activity value of the compartment determined iteratively is above a predetermined or adaptive threshold value.

7. The method of claim 6, wherein compartments having an estimated activity value below a predetermined or adaptive threshold value are excluded from the further iterative procedure.

8. The method of one of claims 1 to 7, wherein step g) includes interrupting the iterative procedure if any of the following criteria are met if:

• the number of active compartments has reached a predetermined value,

• the number of active compartments has not changed during a predetermined number of previous iterative steps,

• a predetermined maximum number of iterations has been reached, and/or

• a predetermined number of local maxima of the estimated activity value has been found in the 3-dimensional volume.

9. The method of one of claims 1 to 8, comprising the step of iteratively reducing the size of the compartments.

10. The method of one of claims 2 to 9 comprising obtaining time-resolved data of the activity of the active sources.

11. The method of one of claims 1 to 10, wherein the physical quantities measured in step a) are EEG and/or MEG data of the human or animal brain, wherein the n detectors are attached to the human scalp or placed around the human or animal head, respectively.

12. The method of claim 11 , wherein the source activities are calculated according to the formula s(t) = T^'1 d(t), wherein s(t) are the time dependent source activities, T^"1 the inverted transfer matrix of the head derived from a volume conductor model of the head, and d(t) are the time dependent measurement results obtained from the n detectors.

13. Use of the method of claim 11 or 12 for the diagnosis of neurolocigal and psychiatric disorders, for example epileptic disorders, stroke, migraine, dementia, depression etc., or for the monitoring of brain function during intensive care or anaestesia.

14. The method of one of claims 1 to 10, wherein step a) comprises measuring detector data from the human or animal body, for example EKG and/or MKG data of the heart.

15. A computer program comprising program code adapted to perform the method of one of claims 1 to 12.

16. Computer program product stored on a computer readable medium comprising program code adapted for performing the method of one of claims 1 to 12.

17. Apparatus for determining a finite number of discrete sources within a 3- dimensional volume comprising:

• a plurality of n detectors for measuring a physical quantity,

• an analog-to-digital converter for converting the detector signals into digital data signals,

• processing means adapted for: a) dividing the 3-dimensional volume into a plurality of discrete compartments, b) selecting a subset of m compartments, c) determining an estimate of an activity value (Q-value) for the m compartments based on the measurement data, d) storing and updating the estimated activity values, e) iteratively repeating steps b), c) and d) with updated subsets of compartments, f) determining a number of compartments as active sources based on the stored and updated estimated activity values.

18. A method of imaging a 3-dimensional volume utilizing a finite plurality of discrete sources, comprising the steps of:

a) dividing the 3-dimensional volume into a finite plurality of compartments, b) providing locations and/or orientations of a defined number of active sources within the 3-dimensional volume, c) measuring a physical quantity at a plurality of n detectors, d) determining the activities of the active sources based on the measurement results, and e) creating a 3-dimensional image of the 3-dimensional volume by calculating the activity for each compartment using a probe source in addition to the active sources.

19. A method of time and frequency dependent pattern analysis of the activities of a plurality of discrete sources within a 3-dimensional volume, comprising the steps of:

a) dividing the 3-dimensional volume into a finite plurality of compartments, b) providing locations and/or orientations of a defined number of active sources within the 3-dimensional volume, c) time-resolved measuring a physical quantity at a plurality of n detectors, d) determining activities of the active sources based on the measurement results, e) deriving parameter vector representations of the three parameters time of occurence, spectral content, and duration of activity of each of the active sources, f) performing a spectral-temporal filtering operation by comparing the parameter vector representations of the measured source activities with a previously stored two-dimensional set of basic signal patterns, or with a previously determined two-dimensional set of signal patterns derived from the complete measurement over time using spectral-temporal decomposition techniques on the source signals, respectively, and g) performing an iterative optimization procedure to determine as analysis result a set of discrete activity events for each active source defined by location, orientation, time of occurence, duration and periodicity or signal frequency, respectively.

20. The method of claim 19, wherein the previously stored set of basic signal patterns consists of a plurality of wavelet functions defined by frequency and duration or the number of cycles at each frequency, respectively.

21. The method of claim 19, wherein the previously stored set of predetermined basic signal patterns consists of a set of specific signal patterns defined by duration and number of cycles of the underlying periodicity pattern.

22. The method of claim 19, 20 or 21 , wherein a discrete weighting function is obtained for each activity event of each active source.

23. Use of the method of one of claims 19 to 22 for the detection of pathological patterns in EEG or MEG data, e.g. slow waves, spikes or seizure discharges in epileptic disorders or stroke, or in EKG or MKG data.