Abstract:
An equivalent dipole localization system comprising, potential measurement means for simultaneously measuring the potentials on a plurality of electrodes attached in a living body, dipole model calculation means for assuming a electric dipole at a given position in the living body having a character of a medium and calculating positions corresponding to the plurality of electrodes formed by the electric dipole, correlation function calculation means for calculating a spacial correlation function of the potential on each electrode from the measurement of the potential measurement means, maximum likelihood estimated error calculation means for calculating an error between a measured value obtained by the potential measurement means and a calculated value obtained by the dipole model calculation means according to calculated values obtained by the correlation function calculation means, and equivalent electric dipole setting means for setting an equivalent electric dipole by obtaining the location, the vector component and the intensity thereof of electric dipole that can minimize the error obtained by the maximum likelihood estimated error calculation means.

Description:
BACKGROUND OF THE INVENTION 
     The present invention relates to equivalent electric dipole localization systems and, more particularly, to an equivalent electric dipole localization system, in which the nerve activities of a living body are represented by an electric dipole for obtaining information about the electric dipole sources and also information about the confidence region or confidence limit of the estimated electric dipole inversely from potential distributions projected on the organic surface. The system localizes equivalent current dipole in living tissue. 
     Heretofore, electroencephalographs, electromyographs, evoked potential adders, etc. have been used to measure potentials evoked on the organic surface by nerve activities of living bodies. Also, recently an equivalent dipole localization method has been proposed, in which activated areas in the living body are estimated by measuring potentials evoked on the organic surface by the nerve activities of the living body. In this method, with stimulation of brain cells in the activated areas, electromotive forces are generated to create potential distributions on the scalp. From such potential distributions, the individual activated areas are corresponded by an electric dipole, and the location and vector component of the dipole are calculated from the potential distributions noted above to estimate the locations of the active brain cells, thereby obtaining numerical representation of the active status of the brain. In a specific proposed method of calculation, a plurality of electrodes are fitted on the organic surface, and potentials generated on the electrodes by the nerve activities of the living body are measured simultaneously. Then, an electric dipole is assumed in a predetermined location in the living body has a medium having a certain character, and the potentials at the locations of the electrodes created by the electric dipole are calculated. Further, the square error between the measured and calculated values obtained for each electrode is obtained, and the location and vector component of the electric dipole that minimize the square error are obtained as those of an equivalent electric dipole. 
     In the above prior art equivalent electric dipole localization method, the error function to be minimized is a least square error function of the measured and calculated potentials obtained with each electrode. Such an error function, however, is effective only when there is no background noise correlation among the individual electrodes. For example, S. Kuriki, M. Murase and F. Takeuchi, &#34;Locating accuracy of a current source of neuromagnetic response: simulation study for a single current in a spherical conductor&#34;, Electroencephalography and Clinical Neurophysiology, 1989, Vol. 73, pp. 499-506, report that errors contained in dipole locations estimated to be greatly electroencephalographically affected amount to 1 cm. Brain waves evoked by visual stimulation may contain, in one waveform, components with close cortical regions such as regions 17 and 18 as sources, and it is very difficult to separate such components and specify the locations thereof with a dipole expressing the individual sources. 
     What are known about the spacial characters of the background brain waves are surprisingly few, but there are known characters in the cerebrum, and similar characters could be estimated about signals measured on the scalp. For example, J. C. De Munck, P. C. M. Vijn and H. Spekreijise, &#34;Random dipoles as a model for spontaneous EEG- and MEG activity&#34;, Advances in Biomagnetism (edited by S. J. Williamson et al., Plenum Press, New York, 1989), show that the spacial characters of background brain waves can be described in terms of a linear relation between variance and inter-electrode distance. If such a spacial correlation is utilized, it will be possible to select better functions as the error function noted above, and the conventional least square error function may be utilized after removing such spacial correlation by filtering. 
     SUMMARY OF THE INVENTION 
     The present invention seeks to solve the problems noted above, and its primary object is to provide an equivalent electric dipole localization system, which permits the location and vector component of an equivalent electric dipole to be obtained with high accuracy by providing, as an error function utilized in the equivalent electric dipole localization method, a new error function based on spacial characters of background brain waves. 
     Another object of the present invention is to provide a method of determining mainly the number of equivalent electric dipoles by providing the confidence region and confidence limit of estimated equivalent electric dipole parameters on the basis of the new error function noted above. 
     According to the present invention, there is provided an equivalent dipole localization system comprising: potential measurement means for simultaneously measuring the potentials on a plurality of electrodes attached in a living body; dipole model calculation means for assuming a electric dipole at a given position in the living body having a character of a medium and calculating positions corresponding to the plurality of electrodes formed by the electric dipole; correlation function calculation means for calculating a spacial correlation function of the potential on each electrode from the measurement of the potential measurement means; maximum likelihood estimated error calculation means for calculating an error between a measured value obtained by the potential measurement means and a calculated value obtained by the dipole model calculation means according to calculated values obtained by the correlation function calculation means; and equivalent electric dipole setting means for setting an equivalent electric dipole by obtaining the location, the vector component and the intensity thereof of electric dipole that can minimize the error obtained by the maximal estimated error calculation means. The equivalent dipole localization system further comprises confidence region or limit calculation means for calculating a confidence region or limit of the equivalent electric dipole obtained by the equivalent electric dipole setting means. 
     Other objects and features will be clarified from the following description with reference to the attached drawings. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a schematic representation of an embodiment of the equivalent electric dipole localization system according to the present invention; 
     FIG. 2 is a flow chart for explaining an example of operation of the system shown in FIG. 1; and 
     FIG. 3 is a view showing an example of output results obtained with the equivalent electric dipole localization system according to the present invention, 
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     The present invention will now be described in detail with reference to the drawings. FIG. 1 is a schematic representation of an embodiment of the equivalent electric dipole system according to the present invention applied to the brain in a head as living body 1. 
     In the use of this system, around 30 electrodes 2, for instance, are attached on the organic surface of a portion under measurement, for instance the head, and potentials evoked on the basis of nerve activities in the brain are measured with potential measurement means 7. Measured potential data from the group of electrodes 2 is coupled through an amplifier 3 and a multiplexer 4 to an analog-to-digital converter (A/D) 5, from which digital measured potential data is coupled through an input port 11 to a computer 6. The computer 6 includes a control unit 6a and a first to a fifth operational unit 6b1 to 6b5. An address and a data bus 8a and 8b are connected to a ROM 9, a RAM 10 and an input and an output port 11 and 12. The ROM 9 and RAM 10 serve as memory means for storing necessary programs for signal processing and also data from the potential measurement means 7. The first to fifth operational units 6b1 to 6b5 of the computer 6 each include data extraction means for extracting evoked potential data and background brain wave data, correlation function calculation means, maximum likelihood estimated error calculation means, equivalent electric dipole setting means, and confidence region/confidence limit calculation means. To the input port 11 is connected an external memory unit 13, in which programs for obtaining equivalent electric dipole are stored, and to the output port 12 are connected a CRT or like display means 14 for displaying results of calculations performed in the computer 6 and a printer 15 for storing data and waveforms displayed on the display means 14. 
     The operation of the embodiment having the above construction will now be described with reference to the flow chart of FIG. 2. 
     Referring to FIG. 2, in a first step ST1 the embodiment of the equivalent electric dipole localization system is initialized by turning on a power source (not shown). 
     In a subsequent second step ST2, brain nerve activity potentials are measured by attaching around 30 electrodes on the head as the living body 1 as shown in FIG. 2. The nerve activity potentials may be those evoked in response to various stimuli, such as electric, optical, acoustical ones, or those in the absence of any stimulus. The measured potential data is coupled through the amplifier 3 and the multiplexer 4 to the A/D 5 and thence as digital data through the input port 11 to the computer 6 to be stored in the RAM 10. 
     In a third step ST3, the first operational unit 6b1 extracts evoked potential data R and background brain wave data η from the measured potential data stored in the RAM 10. Where visually evoked potentials are measured, for instance, samples are collected for a predetermined period of time as a sweep for each electrode whenever a trigger signal is generated from a visual stimulator. These data are used to extract the evoked potential data by taking the arithmetic mean of the individual sweeps, and also they are used directly as the background brain wave data. Denoting the number of electrodes by I and the number of samples in each sweep time by J, the evoked potential data is obtained as I×J matrix r 
     
         r.sub.ij (i=1, . . . , I, J=1, . . . , J 
    
     The background brain wave data may be obtained by subtracting the odd sweep data from the even sweep data. Alternatively, depending on the character of the visual stimulus what is obtained by filtering the measured potential data through a digital filter (such as a low-pass filter, a high-pass filter or a band-pass filter) may be used as the background brain wave data. 
     In a fourth step ST4, the second operational unit 6b2 calculates the spacial correlation functions among all the electrodes with respect to the background brain wave data obtained in the step ST3. The spacial correlation functions may be calculated by using, for instance, formula (1) ##EQU1## where X and Y are measured potentials on electrodes x and y, N is the number of time samples, and S x  and S y  are standard deviations of X and Y. With the number I of electrodes, an I-th order square matrix Q can be obtained. The matrix Q is a symmetric matrix, and if it is a non-negative, it is capable of eigenvalue decomposition as 
     
         Q=  .sup.2  .sup.T                                         (2) 
    
     where   is an orthogonal matrix,   is a diagonal matrix having the eigenvalue of the matrix Q as its component, and T is the transposition of matrix. In a fifth step ST5, a matrix C defined by formula (3) 
     
         C= .sup.-1  .sup.T                                         (3) 
    
     is calculated by utilizing the above eigenvalue decomposition. 
     Then, in the same step ST5 the third operational unit 6b3 calculates potentials r -  generated on the scalp from nerve activity potential sources assumed as an electric dipole. The potentials r -  can be expressed by formula (4) 
     
         r.sup.31.sub.ij =Σφ.sub.j (p.sub.l)S.sub.lj      (4) 
    
     
         (i=1, . . . , I, j=1, . . . , J, l=1, . . . , L) 
    
     where S lj  is the lapsed time on the amplitude an l-th potential source, φ 1  (p l ) is the potential generated on an electrode i by the l-th potential source and has a potential source parameter p l . This function can be written as 
     
         φ.sub.il =Σμ.sub.l,λ φ.sub.i,λ (q.sub.l)(5 
    
     Here, λ=1 represents a radial component of the potential source, and λ=2 and λ=3 both represent tangential components of the potential source. Further, φ i ,λ (q l ) specifies the relationship among the potential source, medium character and potential and has a potential source location parameter q l . 
     In a subsequent sixth step ST6, the fourth operational unit 6b4 calculates an error function between the extracted evoked potential data obtained in the third step ST3 and calculated evoked potential data obtained in the fifth step ST5. This error function is defined by formula (6) given hereinunder. The grounds for utilizing this formula will now be described. 
     Denoting the extracted evoked potential data obtained in the third step ST3 by r ij  and the calculated evoked potential data obtained in the fifth step ST5 by r -   ij , the relation between these two data can be given by formula (6.1) 
     
         r.sub.ij =r.sup.-.sub.ij (p)+η.sub.ij                  (6.1) 
    
     
         (i=1, . . . , I, j=1, . . . , J) 
    
     where p is a parameter of the dipole to be estimated, and η ij  is additive noise and represents background brain wave. If η ij  satisfies 
     
         E(η.sub.ij)=0, E(η.sub.ij, η.sub.i&#39;j&#39;)=Q.sub.iji&#39;j&#39;(6.2) 
    
     the dipole model parameter p which minimizes the error function given as formula (6.3) is the optimal parameter in the sense that the parameter has minimum variance. This is commonly termed Gauss-Marcov theorem. 
     
         H(p)=Σ(r.sub.ij -r.sup.-.sub.ij)Q.sub.iji&#39;j&#39;.sup.-1 (r.sub.i&#39;j&#39; -r.sup.-.sub.i&#39;j&#39;)                                        (6.3) 
    
     
         (i,i&#39;=1, . . . , I; j,j&#39;=1, . . . , J) 
    
     When it is assumed that the background brain waves are the stationary process, the covariance matrix Q iji&#39;j&#39;  of the background brain waves can be estimated using actual brain wave data. The estimation adopts cross correlation function if j≠j&#39; and i≠i&#39;, autocorrelation function if j≠j&#39; and i=i&#39;, spacial correlation function if j=j&#39; and i≠i&#39;, and identically unity if j=j&#39; and i=i&#39;. A case will now be considered, in which the covariance of the background brain waves is given by only the spacial correlation function. In this case, the matrix Q is an first order symmetrical matrix and capable of eigenvalue decomposition. 
     
         Q=  .sup.2  .sup.T                                         (6.4) 
    
     where   is an orthogonal matrix,   is a diagonal matrix given as 
     
          =diag (√λ.sub.1, . . . , √λ.sub.I) 
    
     Here, λ k  (k=1, . . . , I) is the eigenvalue of the matrix Q. Thus, the error function of the formula (6.3) is given as formula (6). 
     
         H(p)=T.sub.r [(CR-CR.sup.- (p)).sup.T (CR-CR.sup.- (p))]   (6) 
    
     where ##EQU2## The matrix C is obtained in the fourth step ST4 (formula (3)). 
     In a subsequent seventh step ST7, the location and vector component of electric dipole that make the error function given by the formula (6) in the sixth step ST6 minimum are obtained (which is a commonly termed inverse problem), and a check is done by the computer 6 as to whether the value of the error function is no higher than the reference value. 
     Here, setting 
     
         Rc≡CR,                                               (7.1) 
    
     and 
     
         Ψ.sub.c ≡CΨ                                  (7.2) 
    
     where [Ψ] il  =Ψ i  (P l ) and Ψ i  (p l ) is given by the formula (5), leads to quite an equivalent problem to the inverse problem algorithm proposed by J. C. De Munck in &#34;The estimation of time varying dipoles on the basis of evoked potentials&#34;, Electroencephalography Clinical and Neurophysiology, Vol. 77, 1990, pp. 156-160. In the inverse problem algorithm proposed by J. C. De Munck, two different linear parameters contained in a dipole model (time function of potential source amplitude and direction of the dipole) and one non-linear parameter (location of the dipole) are handled dividedly, and thus the algorithm comprises three steps. A first step is to optimize the time function, a second step is to optimize the direction parameter, and a third step is to optimize the location parameter, these steps being executed in the mentioned order. 
     If the error function is no greater than the reference value, an eighth step ST8 is executed, in which the location of the electric dipole is shifted by the Marquardt&#39;s method, and then the routine is returned to the fifth step ST5. This operation is iterated until the value of the error function is converged. These operations are all executed in the fourth operational unit 6b4. The Marquardt&#39;s method is one of the non-linear optimization methods and executes iterated calculation to obtain approximation of the solution. In the iterated calculation, it is necessary to make a variable shift for reducing the value of the error function. As an algorithm concerning the direction and extent of the shift, the interpolation-extrapolation method may be utilized. The Marquardt&#39;s method and the interpolation-extrapolation method are described in, for instance, Y. Bard, &#34;Nonlinear Parameter Estimation&#34;, Academic Press Inc., 1974. 
     When the value of the error function is converged to be less than the reference value, thus yielding &#34;YES&#34;, a ninth step ST9 is executed, in which the fifth operational unit 6b5 calculates the confidence region of the estimated electric dipole parameter. Especially, the confidence region of the non-linear parameter q* is by an (N×L)-dimensional ellipsoid Gp given by formula (8). 
     
         Gp=[q:(q-q*).sup.T B(q-q*)≦r.sup.2 ]                (8) 
    
     where the matrix B is an NL-th square matrix, and the individual components are obtainable as 
     
         B.sub.kl =ΣA.sub.ijk (q*)A.sub.ijl (q*)              (8.1) 
    
     
         A.sub.ijk ≡∂r.sup.-.sub.Cij /∂q.sub.k 
    
     where r -   Cij  ≡[CR -  ] ij , r 2  is the value of p % point in χ 2  distribution with a degree of freedom of NL, 
     
         P[χ.sup.2.sub.NL ≦r.sup.2 ]=p(%) 
    
     and N is the number of non-linear parameters contained in an electric dipole. The formula (8) means that proper non-linear parameter is present with a probability of p % inside the ellipsoid given by it. 
     Further, by setting the eigenvalue decomposition of the matrix B defined by the formula (8.1) as 
     
         B=Σα.sup.2.sub.k v.sub.k v.sub.k.sup.T         (8.3) 
    
     the confidence region Gp is 
     
         Gp=[q:Σβ.sup.2.sub.k [(q-q*).sup.T v.sub.k ].sup.2 /r.sup.2 ≦1]                                                (8.4) 
    
     where β 2   k  and v k  are eigenvalue and characteristic vector of the matrix B. 
     The p % confidence limit for the true parameter q i  (i=1, . . . , NL) is given as 
     
         q*.sub.i -δ.sub.i ≦q.sub.i ≦q*.sub.i +δ.sub.i(9.1) 
    
     where δ i  is the maximum value of |q* i  -q -   i  | and given as 
     
         δ.sub.i =[Σ(v.sub.ik /β.sub.k).sup.2 ].sup.1/2(9.2) 
    
     with vi=(v il . . . v i ,NL) T   
     In a subsequent tenth step ST10, in case where there are two or more assumed electric dipoles, a check as to whether there is an overlap of electric dipoles is done according to the values of the confidence region and confidence limit calculated in the ninth step. If there is an overlap of electric dipoles, an eleventh step ST11 is executed, in which the number of dipoles is changed, and the routine is returned to the fifth step ST5. This operation is iterated until there is no overlap of dipoles. 
     When there is no overlap of dipoles, thus yielding &#34;NO&#34;, a twelfth step ST12 is executed, in which the location and vector component of the equivalent electric dipole are stored in the RAM 10 or like memory. Then, in a thirteenth step ST13 the location and vector component of electric dipole stored in the twelfth step ST12 are read out and coupled to the display means 14 and printer 15. 
     FIG. 3 shows output results of the equivalent dipole localization system according to the invention, which are displayed on the printer. The upper three circles are a rear, a side and a plan view of a head model approximation by a sphere. In each circle, the estimated dipole is shown by a vector. The asterisk represents the location of the dipole, and the bar represents the magnitude. In the central part of the Figure, the upper and lower limits of the estimated dipole parameter are shown. These values are obtained as a result of calculation executed with the formulas (9.1) and (9.2). The location of the dipole is expressed as a parameter set(r 0 , θ 0 , φ 0 ), where r 0  represents the radial direction which is expressed as a ratio in which the head radius is unity. Other parameters are shown in the lower part of the FIG. 3. 
     In the above embodiment, the third to thirteenth steps ST3 to ST13 are all executed by the control unit 6a and first to fifth operational units 6b1 to 6b5 of the computer 6. 
     As has been described in the foregoing, according to the invention it is possible to localize an electric dipole with high accuracy. In addition, it is possible to promote elucidation of information processing in the brain by numerically converting not only information about location of lesion thought to be a potential source (for instance location of epilepsia focus) but information about location having excitability by external stimulation (such as optical, graphical, acoustical and electrical stimulation). 
     Further, according to the invention it is possible to obtain the confidence region and confidence limit of each of a plurality of assumed electric dipoles. If two electric dipoles are assumed such that their confidence regions or limits overlap, the number of electric dipoles should be thought to be one. Thus, according to the invention it is possible to provide effective means for checking the number of electric dipoles, which frequently poses problems in the equivalent electric dipole localization method.