Une étude mathématique de la localisation d'une tumeur cancéreuse dans le cerveau humain

C. R. Acad. Sci. Paris, t. 327, SCrie I, p. 601-606, Analyse num&ique/Numerical Analysis

Une 6tude mathhmatique d’une tumeur can&reuse Habib AMMARI a Centre

(Regu

RCsumC.

appliqubes,

of Mathematics,

le 25 mai

Nous

de la localisation dans le cerveau

humain

a, Gang BAO b

de mathkmatiques

b Department

1998

1998,

accept&

dtkrivons

UMR

University aprk

CNRS

of Florida, &vision

une mkthode

7641,

I? co 1e polytechnique,

Gainesville,

le 11 septembre

numtkiique

pour

localiser

cerveau humain. 0 AcadCmie des Sciences/Elsevier,

Mathematical of epileptic

foci

analysis of the in the living

FL

32611,

91128

Palaiseau,

France

USA

1998)

une tumeur can&reuse Paris

determination human brain

dans le

of locations

Abstract.

A new numerical method for determining locations of epileptic foci in the living human brain is derived. The method is based on a low-frequency asymptotic analysis of Maxwell’s equations. The human brain is modeled as a heterogeneous medium, where the electric permittivity, magnetic permeability, and conductivity may all be functions. A current dipole is used to model the epilepsy. The inverse source problem in this context is to determine the current dipole from boundary measurements of the fields. Our method is shown to have good convergence properties. Uniqueness and stability results for the inverse problem are also established. The method and results are expected to jnd applications particularly in magnetoencephalography. 0 AcadCmie des SciencesElsevier, Paris

A bridged

English

Version

Understanding the human brain, the most complex organized structure known to exist, presents a great challenge to the scientific community. The present work is devoted to the study of an inverse source problem that arises in determining locations of epileptic foci in the living human brain. Determining locations of epileptic foci is of great clinical interest. In fact, at present the only way to cure a patient with epileptic foci permanently is to remove a small part of the brain by surgery or by means of gamma-rays. The data are measured via magnetoencephalography (MEG) -a non-invasive technique for investigating neuronal activity of the living human brain [2]. The time resolution for the method is better than 11 ms and the spatial discrimination is 2-3 mm for sources in the cerebral cortex. Note pdsentie

par Philippe

G. CIARLET.

0764~4442/98/03270601 0 AcadCmie des Sciences/Elsevier, Paris

601

H. Ammari,

G. Bao

It is well known in the biomedical engineering literature that epilepsy is caused by electrical discharges that originate from a small volume, and hence may be approximated by a current dipole. Due to the complicated detailed structure of the human brain, it must be modeled by a heterogeneous medium, i.e., the magnetic permeability ,v, the electric permittivity E, and the electric conductivity 0 may be spatial functions. In MEG, one usually deals with frequencies that are below 100 Hz. It should be pointed out that since p N ~0 (the free space magnetic permeability), E - lo5 ~~ (the free space electric permittivity), and Q N 0.3 f12-l m-l, the wave length is about 65 m, i.e., much longer than the diameter of the head. The low frequency assumption is crucial in our analysis. In this Note, we formulate and solve the inverse problem for determining the spatial current source distributions in a heterogeneous medium from boundary measurements of the fields. It is shown that the location of a dipole can be uniquely determined from the tangential components of the magnetic and electric fields on the boundary at a fixed low frequency. A stability result for the inverse problem is also proved. Based on a low-frequency asymptotic analysis, we then derive a constructive numerical method for solving the inverse problem. Our results and analysis indicate that the method is accurate with good convergence properties, and computationally attractive. Our method is further extended in the last two sections. In Section 3, higher order asymptotic methods are also constructed and analyzed. Finally, the variable conductivity case is studied in Section 4, where geometric optics solutions are used.

1. Formulation

du probl&me

Soit (~1, $2, x3) un syst&me de coordonnkes cartksiennes associk B une base orthogonale (el, e2, ea). Soit Ri un ouvert born6 rkgulier de R3 de bord I’ et 0’ le complkmentaire de ?$ dans R3. Nous notons n la normale sortante 2 r. Nous considkrons la propagation d’ondes Clectromagrktiques dans le milieu diklectrique fli isotrope caractkid par sa permittivitk Clectrique E, sa permCabilitk magnCtique p et sa conductivitk Clectrique 0. Nous supposons que les coefficients E et p sont complexes et vkrifient !J?eE > 0 et !J2ep > 0. Nous supposons de plus qu’elles sont de classe C2@, 0 < QI < 1, dans le milieu di6lectrique et nous notons ~0, p. la permittiviti et la permkabilite du vide II”. La conductivite g de Rd est supposke constante et celle du vide est suppoke nulle. Supposant que l’objet diklectrique Ri contient un dipale electrique, les Cquations de Maxwell s’krivent : rot EC”) = -iw p H(“)

dans R3,

rot H(“) = iw E EC”) + 0 EC”) + S(x - x0) p. lim (&IPAE-&jE’) Ix/-+m

dans R3,

(1)

=O,

oti S(x - x0) est la masse de Dirac au point x0 (le point source) et po (la direction du dipSle Clectrique) est un vecteur constant. Notre objet est de d&terminer la position du point x0 et le vecteur p. B partir de la mesure des composantes tangentielles des champs EC”) et H(“) sur I’. Ce problbme est crucial en magne’toenc.kphalographie.Le cerveau humain est consid& comme un milieu dielectrique hCtCrog&ne de coefficients p - ~0, ,c - lo5 EO et g - 0,3 R-l m-l, oit ~0 et co sont les coefficients du vide. Une tumeur can&reuse peut &tre modCli&e par un dip8le Clectrique c&ant une activitC Clectrique beaucoup plus intense que l’activitk Clectrique dans le cerveau d’un sujet normal (au repos). La longueur d’onde est t&s grande devant la taille du cerveau. C’est cette don& qui 602

Une

etude

mathkmatique

de la localisation

d’une

tumeur

can&reuse

dans

le cerveau

humain

va nous permettre de construire une mkthode numkique performante pour localiser la tumeur. Avant de prkenter cette m&hode numkique, nous donnons un resultat d’existence et d’unicitk pour les Cquations de Maxwell avec un dipale Clectrique. Nous avons : LEMME 1.1. - Le syst2me de Maxwell (I) est bien posk. De plus, la solution (EC”), H(“)) classe C2>” sur (P \ (x0)) U R”.

2. R&ultat

d’unicitC

et premikre

est de

approximation

Soit e une solution classique des Cquations de Maxwell

dans Ri :

M(")(e)=rotLrote-w2ce-iwoe=O P

dans Ri.

En intkgrant par parties, il est facile de voir que A n.e+-f-

Lrote.E(“) ‘I rP

An.

W

Nous avons le resultat de complktude TH~ORBME

(2)

suivant :

2.1. - Soit

nip”) ={

u E (C27a(s2i))3; M(“)(u)

= 0 dans Ri . >

Soit un champ EC”) dans H(rot) une solution des kquations de Maxwell ci l’exte’rieur de fi2i avec la condition de radiation de Silver-Miiller. Si 1‘identiti prtkkdente (2) est vraie pour tout champ e dans N(@),

alors les champs E (WI et J$“) = -i

poskes dans R3.

rot I+‘) WC1

sont solutions des kquations de Maxwell (l),

La dkmonstration de ce thCor&me de complCtude utilise le lemme suivant, que I’on dkmontre & l’aide des formules de reprksentation intCgrale et du r&ultat de continuation unique pour l’tquation de Helmholtz : LEMME 2.1. - Si Zesfonctions u E (L2(s2i))3, & rotationnel dans (L”(Q’))”

1

rot-rotu-w2cu--iwDu=f P uAn=O sur r, rotu u-n=0

A n=O

et f E (L(ni))”

ve’rifient :

dans Ri,

sur I-,

sur l?,

alors elles sont n&es.

Remarquons que ce lemme est faux si on supprime la condition aux limites u. n = 0 sur r. 11 suffit pour cela de prendre f = -(w” E + iwcT)grad $ et u = grad $ oit $ E H1(Qi), constante sur r. Maintenant, ti l’aide d’un argument de continuation unique, nous dCmontrons par l’absurde le rksultat d’unicitC suivant : 603

H. Ammari,

C. Bao

THBORPME 2.2. - Soient x0 et x1 deux points duns Ri et po et p1 deux vecteurs constants. Soit (Ej, Hj) pour j = 0, 1, l’unique solution des kquations de Maxwell : rot E(y) = -iw 1-1 H(y) rot ;w

duns R”

= zw E E(i) + (r E(“) + 6(x - x.) pj

lim’(,&H!j’Lfi-

x

Ixl++m

&E;)

duns R3,

= 0. J

Si Eo A nlr = El A n/r et Ho A nlr = HI A nlr, alors x0 = x1 et po = pl.

Remarquons que la mesure de la composante tangentielle du champ Clectrique ou du champ magnetique suffit. Les champs E” et H” sont solutions des equations de Maxwell a l’exterieur de Ri avec la condition de radiation de Silver-Mtiller. Nous n’avons pas besoin de mesurer a la fois la composante tangentielle du champ Clectrique et celle du champ magnetique. Elles sont likes par <
(ml . e(xp)

=

PO

Ir

H(“)

A

nse+T

‘IrP

Lrote.E(“‘)

A

n+O(wm),

W

02 m est un entier nature1 non nul pour tout champ e duns N(@).

Alors,

1~:~;“) - pal = O(w”) et Ixim’ -x01 = O(w”). L’idCe pour demontrer ce resultat de stabilite est de construire des solutions asymptotiques dans le noyau N(fii) compte tenu du fait que w < 1. 11 est facile de voir que lorsque la conductivite est constante, pour tout vecteur constant nous pouvons construire une fonction dans N(@) dont la limite, lorsque w -+ 0, est ce vecteur. Ceci donne Iphm’ - pal = O(w”). L’estimation Ix~~) - x01 = O(w”) est obtenue par un autre choix de solutions asymptotiques. Maintenant, pour calculer le point xa et le vecteur PO, nous allons utiliser l’identite (2). Un choix simple de fonctions dans le noyau N(ni) donnera une methode numerique performante. L’idte est de construire des solutions asymptotiques a l’equation M (“)(u(“)) = 0 dans Ri en la variable w. Soit p(O) unefonction scalaire harmonique dans Ri de classe C2+ et u(l) une fonction vectorielle verifiant rot k rot u(l) = io grad p(O) dans Ri. Nous avons : MC”) [grad p(O) + w u(i)

1= O(w2) dans Sti.

Ce qui donne po . grad ‘p(O)(xa) =

I$“) Ir

A n . grad q(O) + i

JrP

1 rot u(l) . EC”)

A

n + O(w).

Un choix possible de fonctions harmoniques dans fli est zl, z2, zs! et eY’X oh < E C3 est tel que < . [ = 0. Les trois premieres fonctions harmoniques determinent aproximativement le vecteur po (2

604

Une

etude

mathkmatique

de la localisation

d’une

tumeur

can&ewe

dans

le cerveau

humain

w prks) et la qua&i&me (avec trois valeurs diffkrentes de <) now permet de localiser le point xg ?I w p&s aussi. Plus pr&.zidment, soient p:) et xi:’ dCfinis par : PC’ ’ ej = p!’

sr

H(“)

A n. ej + i

. lj eit~‘XF’ =

H(“)

1 rot uil’ . EC“) h n J rP

pour j = 1,2,3,

A n . ,$ eicj’X + ,i r i rot, ~5’) . EC”) A n s

sr

pour j = 1,2,3,

(3)

oti les & E C3 sont distincts deux B deux et sont tels que 0 . & = 0 pour j = 1,2,3 et ui.l’ et v(l) 3 sont tels que 1 rot - rot uj (l) = ige. J CL et 1 . &t+~* rot - rot vi (1) _- 20 P dans Ri. Alors, d’aprtis le thCorkme de stabilitk, nous avons : I&)

- PO] = O(w) et Ix:)

- x0\ = O(W).

Ceci nous donne une m&hode simple et performante pour localiser le point source x0 et determiner le vecteur po. Nous disposons de la valeur de H(“) A n et EC”) A n en un nombre fini de points sur le bord !Z (les capteurs), et ?I une longueur d’onde X = 27r/w grande par rapport ?I l’ouvert Ri. Nous commenSons par calculer II’:) et vi.l). Nous interpollons ensuite les intkgrales qui interviennent dans (3). Le critkre ti satisfaire (Aombre de capteurs et mCthode d’interpolation utiliske) est que ces intkgrales soient calculkes avec une prkcision au moins de l’ordre de w (- 10d2). 3. Approximations

d’ordre

Pour calculer une meilleure

ClevC approximation

de p. et x0, on doit choisir u(l) tel que

div(ia u(l) + E grad cp(‘)) = 0. Plus gCntralement,

u(j) est tel que : rot 1 rot ,(i) = E &-2) + i. ,(i-1) P + ~u(j-~)) = 0 dans fli. i div(igu(j)

dans fii,

A partir de MC”)

[

grad p(O) + w u(l) + . . . + urn-1 ~(~-l)

I

= ()(w”+l)

dms fii,

nous avons

1(x0) =J'l-H(")An.grad p(o) +wu(1) +...+urn-1 U(m-1) 1 u(l)+...+urn-l II(~) . EC”) A n + O(w”). I

~0. grad v(O) + w u(l) + . . . + wmwl u(~-~) [

605

H. Ammari,

G. Bao

Si m = 2, nous.definissons pf’ (2)

‘ej

PO

=WU

et it’

(l)(X$l))

.

ej

par :

+

H(“)

r

A n. ej + w u(l)

s u!l) 3

+

w u!‘) 3

I

$1

. & eic~.X?

= w pf’)

. I

E(“) A n

. v(l)(xt))

pour j = 1,2,3, ~(“1

+ i

I

A n. [in, eitj.x

+ w y(l)]

Jr +i

1 rot v(l) 3 + w v!‘)3 J rP [

1.E(“)

A n

pourj = 1,2,3,

oti vjl) et vi2) sont les fonctions associeesa cp(O)(x) = eiEj.X. Nous avons les estimations suivantes :

(pg) - pal = O(w2) et IX!) - x0] = O(w2). En conclusion, cette m&ode numerique qui s’avere Ctre tres peu couteuse devrait nous permettre de trouver la position du point source, la direction du dipole et son intensite, a une precision de l’ordre de w2 N 10P4.

4. Conductkite

variable

Lorsque la conductivite cr est variable, nous devrons utiliser des fonctions cp solutions de div IJ grad cp = 0. Un choix possible est parmi les fonctions d’optique gkom&rique, construites par Sylvester et Uhlmann dans [3]. D’apres ces auteurs, nous savons que pour tout [ E C3, { . < = 0 et I,$] >> 1, il existe une solution ‘pc de la forme :

PdX>= &

(1+R,(x)>,

oti Re
=

J r

1 rot u(l) . EC“‘) A n. J rP

H(“) A n. grad ‘pc + i

Les details des calculs, present& dans cette Note, sont donnes dans [l].

Rdfhrences bibliographiques [l] Ammari H., Bao G., An inverse source problem for Maxwell’s equations [2] Hamalillinen M., Hari R., Ilmoniemi R.J., Knuutila J., Lounasmaa O.V., and applications to noinvasive studies of the working human brain, Rev. [3] Sylvester J., Uhlmann G., A global uniqueness theorem for an inverse 153-169.

606

in magnetoencephalography, prepublication. Magnetoencephalography-theory, instrumentation, Mod. Phys. 65 (1993) 413-497. boundary value problem, Ann. Math. 125 (1987)