System and method for the tomography of the primary electric current of the brain and of the heart

• US 7,092,748 B2
• Filed: 07/26/2002
• Issued: 08/15/2006
• Est. Priority Date: 02/18/2000
• Status: Expired due to Fees

• Current Assignee: Centro Nacional de Investigaciones Científicas
• Original Assignee: Centro Nacional de Investigaciones Científicas
• Inventors: Bosch Bayard, Jorge Francisco, Morales Aguilera, Frank, Valdes Sosa, Petro Antonio, Virues Alba, Trinidad, Riera Diaz, Jorge Javier, Aubert Vazquez, Eduardo Francisco, Trujillo Barretto, Nelson Jesus, Soler McCook, Jorge Miguel, Fuentes Montero, Maria Elena
• Primary Examiner(s): IMAM, ALI M
• Application Number:
  US10/206,006
• Publication Number:
  US 20030093004A1
• Time in Patent Office: 1,481 Days
• Field of Search: 600/544, 600/545, 600/409, 600/426, 600/407
• US Class Current: 600/407
• CPC Class Codes: A61B 5/0536 Impedance imaging, e.g. by tomography
• CPC Class Codes: A61B 5/24 Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
• CPC Class Codes: A61B 5/7264 Classification of physiological signals or data, e.g. using neural networks, statistical classifiers, expert systems or fuzzy systems neural networks per se G06N3/00; expert systems per se G06N5/00

Abstract

A three-dimensional map of the probability of brain or heart functional states is obtained based on electric or magnetic signals, or a combination of both measured in the surface of body. From signals, the statistical descriptive parameters are obtained and a map of its distribution is calculated. The map is the inverse solution of the problem based upon: a) the restriction of the solution to structures with high probability of generating electric activity using for this restriction an Anatomical Atlas and b) imposing that the solution belong to a pre-specified functional space. The probability is determined that this map belongs to a test group. The spatial and temporal correlations of the map are modeled as well as their dependence on experimental covariables. The resulting probabilities are coded in a pseudocolor scale and they are superimposed on a Anatomical Atlas for their interactive three-dimensional visualization.

Claims

1. Method for the Tomography of the Primary Electric Current (TPEC) of the brain (TPECc) and of the heart (TPECk) comprising the following steps of:

pre-processing of a plurality of signals recorded from a subject for the elimination of artifacts and of non-physiological frequency components;

conversion of a plurality of coordinates of sensors to a reference system of a respective head or torso by means of the application of a visco-elastic deformation, that puts into correspondence the coordinates of a plurality of external anatomical markers of the subject with those of a structural image stored in a computing unit (CU);

calculation of respective electric Kernel and magnetic Kernel linear operators that predict the o(t) that would be produced in the subject by the presence of the PEC j(t) in any part of the brain or heart;

calculation of the TPEC by means of Bayesian Hierarchical Estimation Procedure that determines the primary electric current j(t) of the brain or heart using structural information obtained from an anatomical atlas, metabolic information obtained from functional images, and an assumption that the solution belongs to a space of Besov space of given smoothness or determined by a norm in a Megadictionary;

calculation of descriptive parameters either on the basis of the observations o(t) directly, or of the j(t);

calculation of a probability that all or some of the DP are similar to those of a given test group; and

performing the step wherein the obtained TPECc is coded by means of a pseudocolor scale and is overlaid on the anatomical atlas, a scale being fitted according to the probability previously calculated, whereby only statistically significant anatomical sites are highlighted.

2. The method of claim 1, wherein vectors of descriptive parameters (DP) are calculated that are statistics of models of stochastic processes for o(t) or j(t), wherein instances, of DP include at least one of:

subsets of o(t) or j(t);

subsets of o(t) or J(t) or μ

(t), or γ

(t);

subsets of averages, variances and covariances or higher order o(t) or J(t) or μ

(t), or γ

(t) or of their representation in the basis of a Megadictionary;

subsets of cross polyspectra of o(t) or J(t) or μ

(t), or γ

(t) of their representation in the basis of a Megadictionary;

subsets of parametric multivariate Autoregressive functions fitted to o(t) or J(t) or μ

(t), or γ

(t) or their representation in the basis of a Megadictionary;

subsets of non parametric multivariate Autoregressive functions fitted to o(t) or J(t) or μ

(t), or γ

(t) of their representation in the basis of a Megadictionary;

correction measures, mutual information, Granger causality among components of o(t) or J(t) or μ

(t), or γ

(t) of their representation in the basis of a Megadictionary; and

regression coefficients of any one of the previously described types of DP with stimuli emitted by the system or with the subject'"'"'s responses.

3. The method of claim 1, wherein the probability is calculated that all or part of the TPEC belongs to a test group by means of at least one of a theoretical correction and means of statistical resampling methods that compensate for the effect of multiple spatial and temporal correlations, that are due both to the large number of measurements, as well as a consequence of either the estimation process or the physiology of the brain or heart, the effect of said correlations being to increase a type error I in the identification of brain or heart states, eliminating additionally the effects of the dependence of the estimated maps with experimental covariables of interest.

4. The method of claim 1, wherein test groups are constituted to evaluate at least of of an intra and an interindividual change of brain or heart function, examples of test groups including at least one of:

groups of subjects defined as normal by criteria external to their TPEC, which will be specified implicitly by means of a sample of DP of TPEC of normal subjects contained in a database the objective of statistical analysis being to determine the probability that the DP of the obtained TPEC is obtained from a healthy subject;

groups of subjects defined as having pathologies predefined by criteria external to the TPEC, which will be specified implicitly by means of samples of DP of TPEC of patient with the pathologies of interest defined in a database, the objective of statistical analysis being to determine the probability that the TPEC is obtained from a subject with the specified pathology;

groups determined by methods of Unsupervised Pattern Recognition on the basis of the DP of the TPEC of a sample having in common physiological characteristics of interest, the objective of statistical analysis is to determine the probability that the TPEC is obtained from a subject with said physiological characteristics; and

groups consisting of a sample of DP of TPEC obtained in previous moments from the subject under study, where all of the DPs for each analysis moment constitutes a vector time series and the objective of the statistical analysis of this vector time series has two purpose;

to determine changes with time in the functioning of predetermined parts of the brain or heart as reflected in changes of the DP, as well as to carry out predictions about the probability that a vector of DP will attain certain values (reflecting the functional state of the brain or heart) in a moment posterior to the current examination.

5. The method of claim 1, wherein the following statistical operations are carried out:

Calculation of a mean μ

_x(v), and variance, σ

_X(v), of each Test Group for each one of the parameters and each point of a lattice of the generating volume, said means and variances are be expressed by means of a non parametric regression as a function of control variables;

Transformation towards gaussianity x=T(DP), using a non-linear function to obtain a vector x;

Calculation of the z transform z(x_i)=(x_I−

μ

_x(v))/σ

_X(v);

Calculation of the probabilities that z_max and z_min, the maximum and minimum value of z in a region of interest C respectively, exceeds a threshold U;

Calculation of the following transformation;

6. The method defined in claim 5, wherein the calculation of an element I_B→

A/C that quantifies the directed influence of brain or heart region B on region A that is not due to the activity coming from region C comprises the steps of;

a) defining regions A, B, C as sets of arbitrary voxels r₁ε

_g chosen from the generating volume

b) defining the vector valued time series;

a(t) as the values of γ

_m(t) for all voxels belonging to the set A;

b(t) as the values of γ

_m(t) for all voxels belonging to the set B; and

c(t) as the values of γ

_m(t) for all voxels belonging to the set C

c) defining the matrices;

A_t=[a(t),a(t−

1), . . . , a(t−

k), . . . , a(t−

p)]

B_t=[b(t),b(t−

1), . . . , b(t−

k), . . . , b(t−

p)]

C_t=[c(t),c(t−

1), . . . , c(t−

k), . . . , c(t−

p)]

d) defining the matrices;

7. A Method of obtaining a Tomography of the Primary Electric Current (PEC) from the signals originated either from the brain or the heart using an array of external electrical and biomagnetic sensors and a set of anatomical and medical images, the method comprising the steps of:

a) measuring at least one of the electrical and magnetic signals from a sensor array, defining observations for a common time t, designating an EEG/EKG as o₁(t) and a MEG/MKG as o₂(t) respectively, obtaining as a result a vector time series defined at time instants t_i, I=1, . . . , N_t^I which specify a lattice of time instants for electrophysiological signals ℑ

₁;

b) specifying lattices of sensor positions as the set of Cartesian coordinates that define the position and orientation of each sensor in a predefined common body reference system, said lattices of sensor positions being designated ₁ and ₂ for the EEG/EKG and MEG/MKG respectively;

c) selecting an anatomical image and specifying the Cartesian coordinates, in the common body reference system, of its constituent voxels as the lattice of the volume conductor _v having points r_v;

d) specifying for each point r_v of _v, the probability p(s,r_v) of that voxel belonging to a given tissue type s thus defining an Anatomical Atlas;

e) Labeling each point r_v of the lattice of a volume conductor with a tissue type that has the highest probability, that is s_v=arg max_v(p(s,r_v));

f) specifying for each point r_v, of the lattice of the volume conductor _v, a conductivity value corresponding to the tissue label s_v, said values being taken from a predefined set of conductivities, designated as the conductivity profile σ

;

g) selecting a label b that indicates excitable tissue in order to define;

a probability p(b,r_v) of excitable tissue for each points r_v,

a set _g of N_g points r_g for which said probability is non-zero, and

a probability function p_g=p(b,r_g) defined on _g;

h) calculating an electric K₁ and magnetic K₂ lead field matrices by a vector boundary element method such that, given σ

, ₁, ₂, _g, and _v, the following equations hold;

o₁(t)=K₁f₁(t) and o₂(t)=K₂f₂(t), where f₁(t) and f₂(t) are used interchangeably to denote the PEC;

i) selecting one or more of images selected from the group consisting of the following functional modalities;

functional Magnetic Resonance (fMRI), Positron Emission Tomography (PET), and Single Photon Emission Computed Tomography (SPECT), and defining o₃(t), o₄(t), and o₅(t) observations, respectively;

j) specifying the lattices of voxel positions as the set of Cartesian coordinates, in the common body reference system, for each selected functional image, said lattices being designated ₃, ₄, and ₅ for fMRI, PET, and SPECT respectively;

k) specifying the lattices of time instants t_I, I=1, . . . , N_t^m for which each selected functional image m has been sampled, in a common time scale with ℑ

₁, said lattices being designated ℑ

₃, ℑ

₄, and ℑ

₅ for fMRI, PET, and SPECT respectively;

l) calculating an aggregation operator K_m that expresses our knowledge that the ideal functional indicator f_m(t) associated to a m-th image modality (m=3,4,5) defined on _v and ℑ

_m is modified by the image formation process according to a transformation o_m(t)=(K_m*f_m)(t) which;

reduces a spatial resolution of f_m(t) from that of _v to that of _m, and

reduces the temporal resolution of f_m(t) from that of ℑ

₁, to that of ℑ

_m;

m) calculating the physiological operator H_m that expresses our knowledge that the neural or cardiac tissue activation a(t) defined on _g and _m produces physiological processes associated to the m-th image modality according to the transformation f_m(t)=(H_m·

a)(t);

n) calculating estimates of the PEC by joint estimation of all the f_m(t) and of a(t) for all observed m, said estimation comprising steps of;

i. assigning arbitrary initial values to the f_m(t) and a(t),

ii. iteratively modifying the values of f_m(t) and a(t),

iii. calculating a probability measure p which increases when values of f_m(t) and a(t) simultaneously

Reconstruct the o_m(t) with small error as quantified by a risk function, and

Have a small norm in a given Besov space, and

iv. Continuing step (ii) until the probability measure p does not increase;

o) defining γ

_m(t)=f_m(t), m=1, 2, 3, 4, 5, and γ

₆(t)=a(t);

p) calculating from γ

_m(t) the following quantities;

a vector s^m={s_j^m} that quantifies the magnitude of γ

_m(t) for r_iε

_g, and

a value I<sub>B→

A/C</sub> that quantifies the directed influence of brain or heart region B on region A that is not due to the activity coming from region C;

q) calculating probabilities p_G(s_i^m) and p_G(s<sub>i→

j</sub>^m) that the magnitudes and influences obtained for vector s^m={s_j^m} and I<sub>B→

A/C</sub> are typical of a given reference group G; and

r) coding the values of p_G(s_i^m) and p_G(s<sub>i→

j</sub>^m) by means of a color scale and displaying said codes overlaid on the selected anatomical image highlighting statistically significant regions by thresholding to zero those values beneath a chosen significance level to obtain a statistical parametric map for the Tomography of PEC.

8. The method defined in claim 7, wherein the common reference system is a Talairach coordinate system.

9. The method defined in claim 7, wherein the anatomical atlas is selected from the group consisting of following variants:

a) as an Individual Anatomical Atlas by the extraction of tissue probabilities of anatomical images obtained for the specific subject under study being at least one of the following types;

Computed Axial Tomography (CAT), Magnetic Resonance Images (MRI), and post-mortem sections of the head; and

b) as a Probabilistic Anatomical Atlas by one of a rigid or elastic transformation of the subject'"'"'s volume conductor lattice in order to ensure the best correspondence possible to a canonical volume conductor lattice for which a tissue probability distribution has obtained from a sample of normal or pathological anatomical images.

10. The method defined in claim 7, wherein the computation of the Electric and Magnetic Lead Fields matrices by a vector boundary element method that comprises the steps of:

a) specifying the positions of the sensors for the lattice of EEG/EKG and MEG/MKG modality ₁ and ₂;

b) specifying the positions r_g and r_v in the lattice of the generating volume _g, and of the volume conductor _v respectively;

c) specifying the conductivity profile using an approximation such that each set of tissue labels constitutes a collection of N embedded regions for which a conductivity value is constant, σ

={σ

₁, . . . , σ

_N};

d) calculating numerically values for ohmic current densities j_kV=j_k(r_v) and j_kV^ω

=j_k(r_v,ω

) on each r_v of the lattice of volume conductor _v that belongs to the surfaces limiting the regions by means of the following Linear Algebraic Systems of Equations;

Dc=c_∝

−

Γ

c;

wherein c=(c₁;

. . . ;

C_N);

said c_k={j_vK,V=1, . . . , N_k,k+1} for the EEG/EKG or c_k={j_kV^ω

, V=1, . . . , N_k,k+1} for the MEG/MKG;

N_k,k+1 is a number point r_v that belong to a surface separating the regions k and k+1, wherein matrices Γ and

D are being defined as;

11. The method defined in claim 7, wherein the calculation of the aggregation operator K_m associated to the m-th image modality (i=3, 4, 5) is defined by the following operation

12. The method defined in claim 7, wherein the estimation of the PEC f₁(t) for a given subject comprises steps of:

a) assigning initial values to the following sets of variables;

1) tissue activation a(t)=[a(r_g,t)]<sub>rgε

g</sub> for each point r_g of _g

2) functional indicators f_m(t)=[f_m(r_g,t)]<sub>rgε

g</sub> for each point r_g of _g where f(t)=[f_m(t)]<sub>1≦

m≦

M</sub> denotes the totality of measurement modalities obtained for the subject,

3) hyperparameter of the activation Θ

_a,

4) hyperparameter of the intrinsic noise of the functional indicators Θ

_β

=[Θ

<sub>β

m</sub>]<sub>1≦

m≦

M</sub>, and

5) Hyperparameter of an instrumental noise Θ

_e=[Θ

_em]<sub>1≦

m≦

M</sub>;

b) computing a posteriori probability;

13. The method in claim 12, wherein electrophysiological observations o_m(t) m=1,2 for the estimation of f₁(t) is selected from the group consisting of:

EEG alone, MEG alone, and EEG and MEG together.

14. The method in claim 12, wherein the electrophysiological observations o_m(t) m=3, 4, 5 for the estimation of f₁(t) are selected from the group consisting of:

fMRI alone, PET alone, SPECT alone, fMRI and SPECT, fMRI and SPECT, and PET and SPECT.

15. The method in claim 12, wherein the calculation of the likelihoods I_om(o_m(t)−

(K_m*f_m)(t)|Θ

_em comprise the following steps;

a) selecting the use of either e insensitive polynomial risk or a generalized polynomial risk R_x[x|Θ

_x] for the specific purpose of evaluating non gaussian likelihoods;

b) calculating a discrepancy vector Δ

=o_m(t)−

(K_m*f_m)(t);

c) resetting a(t) as a function of the hyperparameters Θ

_em;

d) calculating R_x^C[Δ

|Θ

_x], the risk of Δ

; and

e) calculating I_om(Δ

|Θ

_em)=exp(−

BR_x^C[Δ

|Θ

_x]) where is a preselected constant.

16. The method in claim 12, wherein the calculation of the a priori probability π

_fm(f_m(t)−

(H_m·

a)(t)|Θ

<sub>β

m</sub> for the functional indicators comprises the following steps;

a) selecting indices m, n, s of the Besov space B_n,s^m, with the purpose of defining the smoothness of the estimated functional indicator, a particular but not exclusive choice being m=1, n=1, s=1 which defines the “

algebra of bumps;

”

b) calculating the discrepancy vector Δ

=f_m(t)−

(H_m·

a)(t);

c) resealing the Δ

as a function of the hyperparameters Θ

<sub>β

m</sub>;

d) calculating ∥

Δ

∥

_Bn,sm, the norm of the Δ

in B_n,s^m, said calculation being carried out by first expanding Δ

as

17. The method in claim 12 wherein the calculation of the a priori probability of the activation π

_a(a(t)|Θ

_a) comprises the following steps;

a) calculating a non decreasing function h(·

) of a probabilistic atlas h(p_g);

b) calculating the norm ∥

a(t)∥

_Bn,sm; and

c) calculating π

_a(a(t)|Θ

_a)=exp(−

Dh(p(s,r_g))−

E∥

a(t)∥

_Bn,sm), where D and E are prespecified constants.

18. The method in claim 16, wherein the a priori probability for the following specific form:

π

_f1(f₁(t)−

(H₁·

a)(t)|Θ

<sub>β

1</sub>)=π

₁(f₁(t)−

Ma(t)|Θ

₁)π

₂(μ

|Θ

₂)

where π

₁ and π

₂ are a priori probabilities specified as in claim 15,

19. The method in claim 18, wherein the a priori probability π

₂(μ

|Θ

₂) defined in claim 17 is calculated by an expression selected from the group of expressions consisting

π

₂(μ

|Θ

₂)=exp(−

∥

Λ

_s·

μ

·

Λ

_m∥

_Bn,sm) where Λ

_s is a diagonal matrix that specifies the degree of smoothness that will be imposed at each point of _g and Λ

_m=W·

P, where W is a prespecified weighting matrix and P is a diagonal matrix containing p_g for all r_g, and

by means of successive evaluation of marginal distributions of π

₂(μ

_g|Θ

₂)=p_g·

N(0,σ

_A²)+(1−

p_g)·

N(0,σ

_B²) with N(0,σ

²) a univariate gaussian distribution with mean zero and standard deviation σ

, σ

_A², and σ

_B² being constants selected to be large and small respectively with respect to the expected variation of the orientations.

20. The method in claim 12, wherein the calculations are carried out for at least one of o₁(t) and o₂(t) in the frequency domain by specifying:

that the electrical signals have been subjected to the Fourier transform, t now denoting frequency;

I_om(o_m(t)−

(K_m*f_m)(t)|Θ

_em)=N_R^2,p(0,Σ

_e), where N_R^2,p(μ

_x,Σ

_x) denotes the multivariate real gaussian distribution with mean and variance (μ

_x,Σ

_x)

π

_fm(f_m(t)−

(H_m·

a)(t)|Θ

<sub>β

m</sub>)=N_R^2,p(0,0), H_m=I, an identity matrix

π

_a(a(t)|Θ

_a)=N_R^2,p(0,Σ

_a)