Cette méthode est évaluée pour un dipôle dans un modèle sphérique en ajoutant un bruit blanc gaussien d'une amplitude de 5% de la valeur maximale de chaque modalité aux données. L'étape 1 donne des erreurs inférieures à 1 mm sur les coordonnées du dipôle, tandis que l'étape 3 donne comme pire erreur 2 mm avec un écart type de 6 mm sur la coordonnée X du dipôle. Ces erreurs donnent à penser que l'on est dans un cas très favorable, et ne sont pas comparées avec les erreurs obtenues en EEG seule.

Dans le même esprit, Kasai et al. [1999] proposent de procéder en distinguant le cas où seul l'EEG mesure un signal significatif du cas où les deux modalités mesurent un signal. Dans le premier cas, le dipôle est contraint à rester radial, et le problème inverse est résolu en EEG seul. Dans le second cas, une première étape consiste à résoudre le problème inverse en contraignant le dipôle à rester tangentiel en MEG, puis dans une deuxième étape, on fixe les paramètres de position du dipôle et les composantes tangentielles du moment dipolaire, et la composante radiale du dipôle est déterminée par résolution du problème inverse EEG. Cette méthode est évaluée pour 10 dipôles dans un modèle sphérique avec différents niveaux de bruit ajouté aux signaux. Elle montre une amélioration par rapport aux modalités traitées séparément.

Huotilainen et al. [1998] utilisent sur des données auditives les deux modalités à la fois séquentiellement et conjointement. Dans un premier temps, les dipôles sont localisés en MEG seule, puis dans un deuxième temps, la position des dipôles étant fixée, soit la MEG et l'EEG sont utilisés conjointement pour déterminer l'amplitude des composantes tangentielles du moment dipolaire, soit l'EEG seule est utilisée, les dipôles ayant alors des composantes radiales et tangentielles.

Dans Stok et al. [1987], les données MEG et EEG sont utilisées conjointement dans un même critère à minimiser de type moindre carrés pondéré. L'évaluation du couplage dans un modèle sphérique et avec différent niveaux de bruit blanc gaussien ajouté aux signaux n'est pas détaillée précisément, mais ils écrivent que l'utilisation simultanée des deux modalités donne toujours un meilleur résultat, le nombre de voies et donc le volume d'information étant plus important lorsqu'on utilise les modalités conjointement.

Dans Mosher et al. [1993], l'étude de la borne de Cramer-Rao dans le cas de données couplées MEG et EEG dans un modèle sphérique pour des dipôles tangentiels montre une amélioration par rapport aux bornes calculées pour les modalités séparées. Toutefois la borne de Cramer-Rao est une borne inférieure de l'erreur en dessous de laquelle un estimateur ne peut descendre (l'estimateur peut se tromper moins sur un cas particulier, mais pas systématiquement moins). Pour montrer que l'on a des estimateurs qui s'approchent de la borne de Cramer-Rao, une simulation par la méthode de Monte-Carlo est réalisée, mais avec des données de MEG seulement.
Le travail le plus complet sur l'évaluation du couplage des données MEG et EEG dans le cas de dipôles ponctuels est certainement celui de Fuchs et al. [1998b]. Dans cette étude, plusieurs concepts sont introduits pour tenir compte des deux modalités exprimées dans des unités différentes :

* Pour ramener les données dans une unité commune, ils proposent de pondérer chaque voie par ses caractéristiques individuelles de bruit : les données sont divisées par une estimation de la variance du bruit à cette voie.

* Comme les données EEG dépendent des conductivités absolues alors que la MEG ne dépend que des rapports entre les conductivités des différents milieux (voir annexe D), ils proposent d'ajuster au mieux ce facteur multiplicatif.

* Enfin ils proposent une nouvelle technique de régularisation permettant d'éviter que les composantes quasi-radiales des sources soient surestimées en modifiant le critère à minimiser :

$$S(\mathbf{d}) = \mathbf{f}(\mathbf{d})^T W \mathbf{f}(\mathbf{d}) + \mathbf{d}_{\text{ori}}^T C(\mathbf{d}_{\text{pos}})^T C(\mathbf{d}_{\text{pos}}) \mathbf{d}_{\text{ori}}$$

où $\mathbf{f}_i = (B_{mes}^i - B_{mod}^i(\mathbf{d}))$, $B_{mes}^i$ désigne les données MEG/EEG mesurées à la voie $i$ et $B_{mod}^i(\mathbf{d})$ celles calculées à la voie $i$ dans le modèle pour des dipôles de positions et d'orientations $\mathbf{d} = (\mathbf{d}_{pos}, \mathbf{d}_{ori})$, $W$ est une matrice diagonale de pondération par voie, et $C(\mathbf{d}_{pos})$ est la matrice de régularisation qui est diagonale dans l'espace des orientations propres. Si $\mathbf{d}_{ori}^j$ est une orientation propre associée à la valeur propre $\sigma_j^2$ alors la valeur propre de $C(\mathbf{d}_{pos})$ associée est $C_j = \sigma_j / B_{mod}(\mathbf{d}_{pos}, \mathbf{d}_{ori})^T W B_{mes}$.

Le couplage de la MEG et de l'EEG est alors évalué en ajoutant un bruit blanc de variance fixée aux données avec 19 montages MEG/EEG différents et 12 dipôles tangentiels dans un modèle sphérique, montrant que la combinaison des deux modalités permet une amélioration de la localisation à toutes les excentricités. Ils trouvent des erreurs de localisation inversement proportionnelles au rapport signal sur bruit. L'erreur de localisation normalisée par rapport au rapport signal sur bruit puis multipliée par 100 est améliorée de 0 mm à 1 mm par le couplage par rapport à la meilleure des deux modalités. La méthode générale dans un modèle réaliste avec des sources ayant des composantes radiales est alors appliquée sur des données évoquées somatosensorielles mais dans ce cas on ne connaît pas la localisation réelle des sources.

Les concepts qui sont introduits dans les méthodes de couplage adaptés aux sources distribuées s'appuient largement sur la connaissance des positions des sources élémentaires qui forment le modèle de sources. Ainsi, dans Baillet et al. [1999], pour chaque position de source, le poids associé à chaque modalité est calculé de façon à minimiser la redondance d'information entre les données MEG et EEG pour cette position de source. De même dans Hasson [1991] une méthode en deux étapes est proposée : le problème est d'abord résolu en
résolvant le problème en utilisant à la fois la MEG et l'EEG dans un premier temps, puis en altérant les poids associés à chaque position de source en fonction de la solution trouvée pour résoudre à nouveau le problème inverse en MEG seule. Ces concepts ne peuvent pas être étendus au problème de la localisation de simples sources dipolaires.

## 4.2 Objectifs

Aucune de ces études (hormis Baillet et al. [1999] pour les sources distribuées) ne propose d'évaluation de l'apport du couplage par simulation dans le cas d'un modèle réaliste. Or le modèle sphérique en MEG est très particulier, puisque les sources radiales y sont complètement muettes. Toutes ces études utilisent des dipôles tangentiels ou exploitent l'absence de signal MEG en contraignant les dipôles à rester tangentiels dans ce cas, ce qui représente une hypothèse forte et difficilement applicable dans le cas de modèles réalistes puisque dans ce cas les dipôles quasi-radiaux contribuent malgré tout au signal. Par ailleurs, dans ces simulations effectuées sur des modèles sphériques, un bruit blanc gaussien indépendant de même amplitude est ajouté aux données simulées. Or dans la réalité, le bruit qui est présent sur les données couplées MEG/EEG est corrélé spatialement (y compris entre les voies MEG et les voies EEG) d'une part, et aboutit à des rapports signal sur bruit différents selon les modalités MEG ou EEG d'autre part.

Dans ce chapitre, nous aborderons plusieurs méthodes permettant de coupler, dans un même problème inverse, les modalités MEG et EEG. Le concept d'ajustement de la conductivité introduit par Fuchs et al. [1998b] sera étudié, et celui de normalisation des données par leurs caractéristiques individuelles de bruit qui peut être vu comme une simplification de la recherche du maximum de vraisemblance décrit dans Lütkenhöner [1998] sera utilisé. Une autre modification du critère inspirée d'un filtrage de Wiener spatial sera également étudiée. Ensuite, l'apport du couplage sera évalué par des simulations avec un modèle réaliste en ajoutant aux données simulées un bruit corrélé spatialement, issu d'enregistrements réels couplés MEG/EEG. Le cas des sources quasi-radiales et quasi-tangentielles sera étudié. En revanche, la régularisation tendant à minimiser la composante quasi-radiale des sources ne sera pas utilisée pour le couplage, les données EEG devant permettre d'éviter une surestimation de cette composante dans le cas du couplage. L'utilité de cette régularisation n'est pas remise en cause, mais l'évaluation de son apport dépasse le cadre de ce chapitre.
## 4.3 Paramètre d'ajustement de la conductivité

Selon Fuchs et al. [1998b], notre mauvaise connaissance des conductivités in vivo des tissus pourrait faire échouer la localisation des générateurs lors de la combinaison de données électriques et magnétiques. Dans les cas les plus simples (conducteur homogène infini, modèle sphérique à une couche), le potentiel est inversement proportionnel à la conductivité du milieu, tandis que le champ magnétique ne dépend pas de la conductivité. Dans le cas général, si l'on fixe les rapports entre les conductivités des différents milieux, on retrouve la même dépendance à un facteur multiplicatif de l'ensemble des conductivités : le potentiel électrique est divisé par ce facteur tandis que le champ magnétique est inchangé. Lorsque l'on dispose de données couplées MEG et EEG, il est donc possible, si l'on connaît les sources, de calculer ce facteur d'ajustement de la conductivité.

L'approche proposée par Fuchs et al. [1998b] consiste à se placer dans un cas idéal où l'on est assez sûr d'avoir de bons rapports signal sur bruit, et un seul dipôle tangentiel pour expliquer les données. En pratique, il est difficile de se placer dans ce cas, c'est pourquoi nous proposons ici d'évaluer ce facteur d'ajustement en le laissant comme un paramètre libre supplémentaire dans l'algorithme de minimisation qui résout le problème inverse. Le critère à minimiser s'écrit d'une manière générale (voir section 1.2.4.1) :

$$S(\mathbf{d}) = \mathbf{f}^t W \mathbf{f}$$

où W est une matrice symétrique définie positive, et dans notre cas :

$$\mathbf{f} = \begin{pmatrix} B \\ V \end{pmatrix} - \begin{pmatrix} B_{\text{mod}}(\mathbf{d}) \\ \rho V_{\text{mod}}(\mathbf{d}) \end{pmatrix} M$$

où ρ est le paramètre d'ajustement de la conductivité. Lorsque le minimum du critère est atteint, on a :

$$\left\{ \begin{array}{l} \left( \begin{pmatrix} B_{\text{mod}}(\mathbf{d}) \\ \rho V_{\text{mod}}(\mathbf{d}) \end{pmatrix} \right)^t W \left( \begin{pmatrix} B \\ V \end{pmatrix} - \begin{pmatrix} B_{\text{mod}}(\mathbf{d}) \\ \rho V_{\text{mod}}(\mathbf{d}) \end{pmatrix} M \right) = 0 \\[1em] (V_{\text{mod}}(\mathbf{d}) M)^t W (\rho V_{\text{mod}}(\mathbf{d}) M - V) = 0 \end{array} \right.$$

La première équation donne :

$$M(\rho) = (B_{\text{mod}}^t(\mathbf{d}) W B_{\text{mod}}(\mathbf{d}) + \rho^2 V_{\text{mod}}^t(\mathbf{d}) W V_{\text{mod}}(\mathbf{d}))^{-1} (B_{\text{mod}}^t(\mathbf{d}) W B + \rho V_{\text{mod}}^t(\mathbf{d}) W V)$$

et la deuxième donne :

$$\rho = \frac{(V_{\text{mod}}(\mathbf{d}) M(\rho))^t W V}{(V_{\text{mod}}(\mathbf{d}) M(\rho))^t W (V_{\text{mod}}(\mathbf{d}) M(\rho))}$$
Ce qui nous donne une équation pour le paramètre $\rho$ que l'on peut résoudre par dichotomie par exemple. Pour éviter de nombreux calculs d'inverses de matrices dans ce cas, il est utile d'utiliser la matrice inversible $X$ qui diagonalise les deux matrices symétriques positives (cf Golub et Van Loan [1996], p 465) :

$$ egin{cases} B_{mod}^{t}(\mathbf{d})WB_{mod}(\mathbf{d}) = X^{t}S_{B}X, \text{ et} \\ V_{mod}^{t}(\mathbf{d})WV_{mod}(\mathbf{d}) = X^{t}S_{V}X \end{cases} $$

donne :

$$ (B_{mod}^{t}(\mathbf{d})WB_{mod}(\mathbf{d}) + \rho^{2}V_{mod}^{t}(\mathbf{d})WV_{mod}(\mathbf{d}))^{-1} = X^{-1}(S_{B} + \rho^{2}S_{V})^{-1}X^{-t} $$

où la matrice à inverser est diagonale, donc il suffit d'inverser un à un chacun de ses coefficients.

Cette méthode permet de calculer $\rho$ dans la partie “linéaire” de l’algorithme de minimisation. Une autre approche, plus facile à implémenter, consiste à ajouter le facteur $\rho$ comme un paramètre non-linéaire, ce qui présente l’inconvénient d’augmenter le nombre de paramètres à optimiser, donc de ralentir l’algorithme, et d’augmenter le risque de tomber dans un minimum local. C’est néanmoins l’approche que nous avons utilisée pour une première validation de cette technique.

Nous avons testé cette méthode sur des données simulées sur un modèle réaliste avec deux dipôles radiaux ou deux dipôles tangentiels, et $\rho = 2$. L’algorithme de minimisation retrouve exactement $\rho = 2$ dans les deux cas. Toutefois le même modèle est utilisé pour créer les données simulées et le problème inverse. En ajoutant du bruit gaussien non corrélat spatiallement sur les données, avec un rapport signal sur bruit de 20, on obtient des valeurs pour $\rho$ comprises entre 1.9 et 2.1. Bien que ces évaluations soient difficilement généralisables, puisqu’elles ne portent que sur deux configurations de dipôles et un modèle de bruit peu réaliste, elles indiquent que la marge d’erreur sur des données réelles risque d’être importante.

Pour illustrer cette difficulté, nous avons testé cette approche dans le cas de données réelles auditives provenant de Yvert et al. [2001a] (sons de 1000 Hz dans l’oreille droite). Il s’agit d’un enregistrement de deux fois 37 voies MEG (gradiomètre axial BTi placé successivement de manière à couvrir l’hémisphère gauche et l’hémisphère droit) et 32 voies EEG (Neuroscan) dont 26 sont retenues ici. A chaque latence, une solution au problème inverse composée de deux dipôles (un dans chaque hémisphère) a été recherchée, et la valeur du paramètre d’ajustement de la conductivité a été relevée pour cette solution. Pour éviter de tomber dans un minimum local, les paramètres des dipôles ont été initialisés à partir d’une première localisation en MEG seule, puis la position des dipôles étant fixée, l’orientation des dipôles a été recherchée en EEG seul avec $\rho = 1$.

Les valeurs trouvées pour ce paramètre d’ajustement de la conductivité sont représentées figure 4.1. En l’absence d’ajustement de la conductivité, ce paramètre vaut 1. L’erreur
<figure><img src="image_55.png" /><figcaption>FIG. 4.1: Paramètre d'ajustement de la conductivité trouvé pour chaque latence sur des données auditives. 'ses1', 'ses2', et 'ses3' correspondent à trois sessions d'enregistrements effectués sur le même sujet, mais des jours différents.</figcaption></figure>

résiduelle est inférieure à 10% entre 50 et 120 ms. Or, le paramètre d'ajustement de la conductivité varie d'une latence à l'autre et d'une session d'enregistrement à l'autre de manière importante : entre 0.7 et 1.6. Il y a plusieurs explications à cette variabilité : les données ne sont peut-être pas dûes à simplement deux sources dipolaires, le nombre de paramètres pour l'algorithme de minimisation (Marquardt) est important (deux dipôles et $ρ$), ce qui peut dégrader la qualité de l'optimum obtenu, enfin il s'agit de données réelles, donc bruitées, et pouvant présenter des variabilités physiologiques entre sessions. On notera toutefois qu'un grand nombre des intersections entre les courbes des différentes sessions entre 85 et 105 ms se situent entre 1.1 et 1.2.

Pour pouvoir utiliser cette technique, il faut donc au préalable mettre au point un protocole permettant d'obtenir avec fiabilité une valeur de $ρ$. Fuchs et al. [1998b] propose de se placer dans un cas où l'on a un dipôle tangentiel et un bon rapport signal sur bruit. En pratique, cependant, il est difficile de s'assurer que l'on est bien dans ce cas. De plus, dans le cas d'un modèle sphérique, la valeur optimale trouvée pour ce paramètre à partir d'un dipôle situé dans le cortex somesthésique risque d'être affectée par une compensation de l'erreur locale sur l'épaisseur des différentes couches. Cette valeur optimale pourrait ne pas convenir pour des sources temporales où le crâne a une épaisseur différente.
## 4.4 Modèle de bruit

L'un des buts du couplage MEG/EEG est de réduire l'incertitude sur l'estimation des sources qui peut être appréciée par l'intervalle de confiance des solutions trouvées. Pour évaluer cet intervalle de confiance, nous allons ajouter des réalisations de bruit réaliste à des données simulées. Pour chaque réalisation de bruit, l'algorithme de résolution du problème inverse sera appliqué, et les localisations ainsi obtenues permettront de définir un nuage de points dont on pourra apprécier la taille. Pour cela, il nous faut nous donner un modèle de bruit aussi réaliste que possible.

### 4.4.1 Comment générer un bruit réaliste ?

Il est courant d'ajouter un bruit gaussien indépendant à chaque point de mesure. Ce modèle de bruit pose problème car en réalité le bruit est corrélé spatialement. Une possibilité pour créer du bruit plus réaliste présentant une corrélation spatiale est présentée dans De Munck et al. [1992]. Il s'agit de simuler l'activité EEG/MEG spontanée grâce à un certain nombre de sources dipolaires, dont la position, l'orientation et l'amplitude sont régis par des lois stochastiques. Toutefois ce modèle est difficile à paramétrer : choix du nombre de sources, des lois qui régissent position, orientation et amplitude des dipôles.

Pour générer un bruit réaliste, nous sommes partis d'un enregistrement de MEG et d'EEG spontané que nous avons moyenné pour obtenir deux types de bruits (bruit1 et bruit2) présentant des corrélations spatiales différentes (voir figure 4.2). L'obtention de ce bruit réaliste est décrite en détail en annexe I.

<figure><img src="image_56.png" /><figcaption>FIG. 4.2: Représentation graphique de l'intercorrélation normalisée entre les voies pour bruit1 à gauche et bruit2 à droite. Les 29 premières voies sont des voies d'EEG (carré en bas à gauche), les 151 dernières voies sont des voies de MEG (carré en haut à droite). Pour une abscisse i et une ordonnée j, la couleur représente la valeur absolue de l'intercorrélation normalisée $\left| \frac{\sigma_{ij}}{\sigma_i \sigma_j} \right|$ entre la voie i et la voie j. La diagonale ascendante correspond à l'autocorrélation normalisée, qui vaut 1.</figcaption></figure>
## 4.5 Modèle utilisé pour la génération des signaux bruités

Le modèle géométrique utilisé pour cette simulation était un modèle à trois couches composées de 3620 triangles pour le scalp, 2800 pour la surface externe de l'os, et 2240 pour la surface interne. Il est présenté figure 4.3. Pour ce modèle, la méthode des intégrales de surface avec interpolation linéaire a été utilisée. Pour les différentes configurations de sources envisagées, le problème direct a été calculé. À ce signal simulé on a ajouté les réalisations de bruit décrites dans la section précédente, en ajustant la contribution du bruit de façon à obtenir les rapports signal sur bruit souhaités, ce qui est équivalent à faire varier l'amplitude des moments dipolaires.

<figure><img src="image_57.png" /><figcaption>FIG. 4.3: Modèle utilisé pour la méthode des intégrales de surface.</figcaption></figure>

## 4.6 Évaluation de l'apport du couplage

### 4.6.1 Utilisation d'une méthode rapide pour le calcul des problèmes directs

Pour chaque position des sources, il y avait autant de problèmes inverses à résoudre que de réalisations de bruit, ce qui nous a conduit à résoudre au total 197400 problèmes inverses. Pour y parvenir, il nous a fallu une méthode rapide de résolution du problème inverse. Ce besoin de méthode rapide de résolution du problème inverse existe également en pratique, pour la modélisation de données auditives par exemple. C'est pourquoi nous avons développé avec Blaise YVERT une méthode basée sur l'interpolation de champs engendrés précalculés sur une grille. Cette méthode est détaillée et évaluée en annexe J.
## 4.6.2 Montages utilisés

En couplant les données MEG et EEG, on utilise plus de capteurs que pour une seule modalité. Pour s'affranchir de l'effet propre au nombre de capteurs et évaluer réellement l'apport du couplage MEG/EEG lié à la complémentarité des données, nous avons utilisé plusieurs montages avec différents nombres de voies MEG/EEG. ainsi certains montages MEG/EEG avaient moins de capteurs que certains montages correspondant à une seule modalité. Ces différents montages étaient les suivants :

* EEG à 29 électrodes (eeg29)

* EEG à 64 électrodes (eeg64)

* EEG à 151 électrodes (eeg151)

* EEG à 256 électrodes (eeg256)

* MEG avec 151 voies (meg)

* EEG à 29 électrodes et MEG avec 151 voies (megeeg29)

* EEG à 64 électrodes et MEG avec 151 voies (megeeg64)

Nous n'avons utilisé qu'un seul montage de MEG car notre idée était de rechercher, lorsqu'on dispose d'une machine de MEG, si il était utile d'ajouter des électrodes et combien. Par ailleurs une limitation de nos programmes informatiques nous interdisait de mettre plus de 256 voies au total. Les positions des électrodes dans les différents montages sont représentés figure 4.4.

## 4.6.3 Critères à minimiser

La résolution du problème inverse consiste à rechercher le minimum d'un critère qui mesure la distance entre le signal théorique engendré par une configuration de sources et le signal mesuré. Dans l'esprit du maximum de vraisemblance, mais sans prendre en compte la covariance spatiale pour ne pas trop alourdir les calculs, nous avons choisi comme premier critère :

$$ \text{criterel}(\mathbf{d}) = \sum_{i=1}^{\text{Nb voies}} \frac{(V_i - \hat{S}_i(\mathbf{d}))^2}{\sigma_{B_i}^2} $$

où d représente les paramètres (position, orientation, amplitude) des sources dipolaires dont le nombre est supposé connu, $V_i$ représente le signal mesuré à la voie i, $\hat{S}_i(\mathbf{d})$ représente le signal théorique engendré par des sources de paramètres d, et $\sigma_{B_i}^2$ est une estimation de la variance du bruit $B_i$ ($V_i = S_i + B_i$) à la voie i.
<figure><img src="image_58.png" /><figcaption>FIG. 4.4: Positions des électrodes dans les configurations à 29, 64, 151 ou 256 électrodes et des spires collectrices pour la MEG.</figcaption></figure>

Le deuxième critère que nous avons utilisé est inspiré du filtre de Wiener. Nous l’avons introduit dans le but de pondérer chaque voie en fonction du rapport signal sur bruit, afin de limiter l’influence du bruit sur la solution. Ce deuxième critère s’écrit :

$$ \text{criter2}(\mathbf{d}) = \sum_{i=1}^{\text{Nb voies}} \frac{(\alpha_i V_i - \hat{S}_i(\mathbf{d}))^2}{\sigma_{B_i}^2} $$

où $\alpha_i$ est choisi pour minimiser la moyenne $\delta = (\alpha_i V_i - S_i)^2$.

Si l'on fait l'hypothèse d'un bruit additif $B_i$ de moyenne nulle, le signal $S_i$ à la voie i s'écrit $V_i = S_i + B_i$.

On a

$$ \delta = \overline{\alpha_i^2 V_i^2} - 2\alpha_i V_i S_i + S_i^2 $$

et donc

$$ \delta = \alpha_i^2 \overline{V_i^2} - 2\alpha_i \overline{V_i} S_i + S_i^2 $$

Ce polynôme du second degré atteint son minimum quand :

$$ \alpha_i = \frac{S_i \overline{V_i}}{\overline{V_i^2}} $$
et comme $\overline{V_i} = S_i$ et $\overline{V_i^2} = S_i^2 + \overline{B_i^2}$,

$$ \alpha_i = \frac{S_i^2}{S_i^2 + \overline{B_i^2}} $$

Ce que l'on peut également écrire :

$$ \alpha_i = \frac{1}{1 + \frac{B_i^2}{S_i^2}} $$

Notons que $0 < \alpha_i < 1$, et que pour un bruit de moyenne nulle, $\overline{B_i^2} = \sigma_{B_i}^2$.

Pour comprendre intuitivement l'effet de ce coefficient $\alpha_i$, on peut considérer les différents cas de figure suivant pour une voie $i$ donnée :

* Si le signal $S_i$ est faible et si le rapport signal sur bruit est faible, alors $\alpha_i$ est proche de 0. La contribution de la voie $i$ au critère tend alors à minimiser le signal modèle sur cette voie, mais cela a peu d'importance car $S_i$ est petit, et cela évite de "modéliser" le bruit.

* Lorsque le bruit $B_i$ est important et que le rapport signal sur bruit est petit, $\alpha_i$ est proche de 0, et la variance du bruit $\sigma_{B_i}^2$ est grande, donc cette voie compte peu dans le critère, ce qui permet une fois encore d'éviter de "modéliser" le bruit.

* Lorsque le bruit $B_i$ est faible et que le rapport signal sur bruit est grand, alors $\alpha_i$ est proche de 1 et la voie $i$ est alors prise en compte pleinement.

Ces deux derniers cas de figure sont en fait également pris en compte par le critère 1. En revanche, le premier cas de figure est mieux pris en compte par le critère 2 qui ne s'appuie pas uniquement sur la variance du bruit $\sigma_{B_i}^2$ mais sur le rapport signal sur bruit : $S_i^2/\overline{B_i^2}$.

Pour chacun de ces critères, il faut estimer la variance du bruit $\sigma_{B_i}^2$, ce qui peut être fait par exemple à partir de ce qui est enregistré dans la période préstimulus. Dans notre cas, puisque c'est une simulation, nous avons pu utiliser la valeur exacte de cette variance. De plus, pour le deuxième critère, il faut estimer le rapport signal sur bruit.

### 4.6.3.1 Estimation du rapport signal sur bruit

Le rapport signal sur bruit (signal to noise ratio, SNR) à la voie $i$ est tel que :

$$ SNR_i^2 = \frac{S_i^2}{\overline{B_i^2}} $$

$\overline{B_i^2} = \sigma_{B_i}^2$ pour un bruit de moyenne nulle. Comme on ne connaît pas $S_i$, mais $V_i$, on utilise la formule $S_i = V_i - B_i$ pour trouver :

$$ SNR_i^2 = \frac{(V_i - B_i)^2}{\overline{B_i^2}} $$
Ce qui donne, en développant :

$$SNR_i^2 = \frac{\overline{V_i^2} - 2\overline{V_i}\overline{B_i} + \overline{B_i^2}}{\overline{B_i^2}}$$

Et comme $\overline{V_iB_i} = (S_i - B_i)\overline{B_i} = S_i\overline{B_i} - \overline{B_i^2} = -\overline{B_i^2}$,

$$SNR_i^2 = \frac{\overline{V_i^2} - \overline{B_i^2}}{\overline{B_i^2}}$$

et finalement :

$$SNR_i^2 = \frac{\overline{V_i^2}}{\overline{B_i^2}} - 1$$

C'est cette formule que nous utiliserons pour calculer

$$\alpha_i = \frac{1}{1 + \frac{\overline{B_i^2}}{\overline{S_i^2}}}$$

### 4.6.4 Critères de qualité adoptés pour l'évaluation

#### 4.6.4.1 Erreur de localisation

Le premier critère de qualité important pour cette évaluation est l'erreur de localisation. Dans le cas où un seul dipôle est responsable du signal, l'erreur de localisation est la distance entre le dipôle à l'origine du signal, et le dipôle obtenu par résolution du problème inverse.

#### 4.6.4.2 Erreur d'orientation

Un deuxième critère important est l'erreur d'orientation entre l'orientation du dipôle à l'origine du signal et celle du dipôle équivalent obtenu par résolution du problème inverse. En effet, dans deux régions anatomiquement proches, les cellules pyramidales de deux groupes de neurones peuvent être orientées différemment, et l'orientation du dipôle peut ainsi permettre de discriminer deux zones actives.

#### 4.6.4.3 Erreur sur l'amplitude du moment dipolaire

Un troisième critère est l'erreur sur l'amplitude du moment dipolaire. L'amplitude du moment dipolaire donne une indication sur le nombre de neurones impliqués.
<figure><img src="image_59.png" /></figure>

### 4.6.4.4 Taille de l'ellipsoïde de dispersion

Un quatrième critère est la stabilité de la solution trouvée au bruit. Pour chaque position de dipôle, les solutions au problème inverse pour les différentes réalisations de bruit forment un nuage de points. Nous avons défini l'ellipsoïde de dispersion ayant pour centre le centre de gravité du nuage, pour axes les axes principaux du nuage, et pour excentricités les écarts-types de la dispersion du nuage selon chaque axe. Le volume de cette ellipsoïde permet d'apprécier la stabilité de la solution au bruit.

### 4.6.4.5 Biais de localisation

Un cinquième critère est le biais de localisation, c'est à dire la distance entre le centre de l'ellipsoïde de dispersion et le dipôle à l'origine du signal.
## 4.7 Première configuration : sources à plusieurs profondeurs

Dans un premier temps, nous avons voulu évaluer l'apport du couplage en fonction de la profondeur de la source. Nous avons placé une source dipolaire à la fois sur un axe temporal, à 11 profondeurs différentes (voir figure 4.5). Pour chaque profondeur, deux orientations ont été envisagées, l'une selon l'axe x, ce qui donne des sources à peu près radiales, et l'autre selon l'axe y, ce qui donnait des sources à peu près tangentielles. Les réalisations de bruit ont été

<figure><img src="image_60.png" /><figcaption>FIG. 4.5: Première configuration des sources : pour 11 profondeurs sur un axe temporal, les sources étaient orientées soit selon l'axe x, soit selon l'axe y.</figcaption></figure>

multipliées par un même facteur de 0.01 quelle que soit la profondeur de la source, ce qui permettait d'avoir un rapport signal sur bruit qui varie en fonction de la profondeur. Les rapports signal sur bruit ainsi obtenus sont représentés figure 4.6.

### 4.7.1 Résultats

#### 4.7.1.1 Effet du nombre d'électrodes pour le couplage

L'une des questions qui se pose dans le couplage des deux modalités MEG et EEG est le nombre d'électrodes nécessaires. Plus généralement, on peut se demander d'une part si, disposant d'une MEG, il est utile d'enregistrer en même temps un petit nombre de voies EEG, ou si un nombre plus grand ne serait pas nécessaire, et d'autre part si on n'aurait pas obtenu une précision comparable avec un grand nombre de voies EEG sans la MEG.

Les différents critères d'appréciation de l'erreur de localisation pour la première configuration de sources et bruit1 sont représentés figure 4.8.
<figure><img src="image_61.png" /><figcaption>FIG. 4.6: Rapports signal sur bruit obtenus pour la première configuration de sources.</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>FIG. 4.7: 4 réalisations de signal bruité pour une source dipolaire profonde quasi-tangentielle placée sur un axe temporal. En haut : EEG (SNR = 6), en bas : MEG (SNR=7).</figcaption></figure>
<figure><img src="image_63.png" /><figcaption>FIG. 4.8: Pour la première configuration de source (voir figure 4.5) et bruit1, de haut en bas : l'erreur de localisation D(mm), l'erreur d'orientation θ(°), l'erreur d'amplitude du moment dipolaire ΔM(%), la taille de l'ellipsoïde (mm³), et le biais de localisation (mm). Les barres d'erreur correspondent à une erreur standard. Les courbes rouges correspondent à l'EEG, les courbes bleues à la MEG et les courbes bleues au couplage MEG/EEG. Les symboles pleins correspondent à critere1, et les symboles ouverts à critere2.</figcaption></figure>
<figure><img src="image_64.png" /><figcaption>FIG. 4.9: Pour la première configuration de sources (voir figure 4.5) et bruit2, de haut en bas : l'erreur de localisation D(mm), l'erreur d'orientation θ(°), l'erreur d'amplitude du moment dipolaire ΔM(%), la taille de l'ellipsoïde (mm³), et le biais de localisation (mm). Les barres d'erreur correspondent à une erreur standard.</figcaption></figure>
* Dans le cas des dipôles orientés selon y (orientation plutôt tangentielle), et quel que soit le type de couplage envisagé, critère1 ou critère2, avec 29 ou 64 électrodes, le couplage donne de meilleurs résultats que l'une ou l'autre des modalités seules, que ce soit en terme d'erreur de localisation, d'orientation, ou en terme de taille d'ellipsoïde. Pour l'erreur sur le moment et le biais de localisation, les résultats sont comparables à ceux obtenus pour les modalités seules. Un plus grand nombre d'électrodes dans le couplage permet une plus grande précision dans le cas des dipôles profonds, pour les dipôles périphériques le nombre d'électrodes n'influe pas sur la précision.

* Dans le cas des dipôles orientés selon x (orientation plutôt radiale), l'avantage du couplage ne ressort pas aussi clairement. Pour l'erreur de localisation, avec le critère1, le couplage avec 29 électrodes donne une précision un peu moins bonne que l'EEG seul. Toutefois, pour les dipôles profonds, eeg29 donne même de meilleurs résultats qu'eeg64, ce qui laisse supposer un effet du hasard dans ce cas. Le couplage avec 64 électrodes donne une précision comparable à celle de l'EEG pour cette orientation, voir même un peu meilleure dans le cas du critère2.

Dans cette orientation plutôt radiale des sources, les erreurs de localisation obtenues avec la MEG seule sont plus importantes, de même que pour les sources les plus profondes avec l'orientation selon y. Toutefois, il apparaît que :

* dans les cas où la MEG donne de mauvais résultats (sources radiales ou profondes), le couplage donne une erreur de localisation comparable ou meilleure que l'EEG (sauf pour le cas particulier des dipôles profonds radiaux avec eeg29 et critère1),

* dans les cas où la MEG donne de meilleurs résultats que l'EEG, le couplage donne une erreur de localisation comparable ou meilleure que la MEG.

L'erreur d'orientation, l'erreur sur le moment, et la taille de l'ellipsoïde sont également soit comparables soit meilleures que pour la meilleure des deux modalités avec le couplage. Seule la distance entre l'ellipsoïde et le dipôle (biais) dans le cas du critère2 semble faire exception à cette règle. *Un plus grand nombre d'électrodes dans le couplage permet une meilleure précision dans les cas où la MEG ne donne pas de bons résultats, mais l'apport du couplage apparaît même avec un petit nombre d'électrodes.*

La même simulation a été réalisée avec bruit2, les résultats sont représentés figure 4.9. Cette simulation confirme les résultats de la simulation précédente : *le couplage donne une précision au moins égale à la meilleure des deux modalités MEG ou EEG, et ce même avec un petit nombre d'électrodes*. Dans les simulations suivantes, pour simplifier la présentation des résultats, nous ne garderons que les modalités meg, eeg64 et le couplage megeeg64.
### 4.7.1.2 Choix d'un critère pour le couplage

Les figures 4.8 et 4.9 montrent systématiquement une amélioration de l'erreur de localisation, de l'erreur d'orientation et de la taille de l'ellipsoïde avec critère2 par rapport à critère1. En revanche, l'erreur sur l'amplitude du moment dipolaire est détériorée par critère2, et un biais de localisation apparaît. L'étude du comportements des deux critères pour les orientations quasi-radiale et quasi-tangentielles nous permettra de mieux comprendre leurs apports et leurs limites.

## 4.8 Deuxième configuration : sources quasi-radiales et quasi-tangentielles

Pour deux profondeurs sur le même axe temporal, nous avons placé des sources selon les trois orientations propres quasi-radiale et quasi-tangentielles définies section 1.2.2. Les valeurs propres obtenues pour le dipôle le plus profond sont 1.05 pour l'orientation quasi-radiale et 45.9 et 50.9 pour les orientations quasi-tangentielles, et pour le dipôle le plus périphérique, 0.86, 625 et 702. Les positions de ces dipôles sont représentées figure 4.10.

<figure><img src="image_65.png" /><figcaption>FIG. 4.10: Deuxième configuration des sources : 2 profondeurs * 3 orientations (rouge : quasi-radiales, jaune et bleu : quasi-tangentielles) sur un axe temporal.</figcaption></figure>

Les réalisations de bruit ont été multipliées par un facteur différent pour chaque source de façon à avoir pour l'EEG à 29 voies un rapport signal sur bruit de 5 (voir figure 4.11).
<figure><img src="image_66.png" /><figcaption>FIG. 4.11: Rapports signal sur bruit obtenus pour la deuxième configuration de sources. rad1, tan1a et tan1b correspondent à la position de source la plus profonde pour les orientations quasi-radiale et quasi-tangentielles. rad2, tan2a, tan2b correspondent à la position de source la plus périphérique pour les orientations quasi-radiale et quasi-tangentielles.</figcaption></figure>

## 4.8.1 Résultats

### 4.8.1.1 Apport du couplage

Les résultats obtenus avec la deuxième configuration de sources sont représentés figure 4.12. On y constate que le couplage permet d'obtenir une meilleure précision de localisation et d'orientation que l'une ou l'autre des modalités, et l'amplitude du moment dipolaire est estimée avec une précision comparable à la meilleure des deux modalités, quel que soit le critère utilisé. La taille de l'ellipsoïde est également réduite par le couplage, en revanche le biais de localisation peut être légèrement détérioré pour les orientations radiales des sources, en particulier avec critère2.

### 4.8.1.2 Choix d'un critère en MEG seule

Les résultats obtenus avec la deuxième configuration de sources, confirment les observations de la section précédente : dans le cas de la MEG seule et des orientations quasi-radiales des sources, le critère2 permet d'obtenir une meilleure localisation et une meilleure orientation des sources que le critère1, au prix d'un moment dipolaire largement sous-estimé. On constate des erreurs d'orientation plus importantes en MEG avec critère2
<figure><img src="image_67.png" /><figcaption>FIG. 4.12: Pour la deuxième configuration de sources (voir figure 4.10), de haut en bas : l'erreur de localisation D(mm), l'erreur d'orientation θ(°), l'erreur d'amplitude du moment dipolaire ΔM(%), la taille de l'ellipsoïde (mm³), et le biais de localisation (mm). Les barres d'erreur correspondent à une erreur standard. rad1, tan1a et tan1b correspondent à la position de source la plus profonde pour les orientations quasi-radiale et quasi-tangentielles. rad2, tan2a, tan2b correspondent à la position de source la plus périphérique pour les orientations quasi-radiale et quasi-tangentielles.</figcaption></figure>
pour les orientations tangentielles des sources périphériques (tan2a et tan2b) par rapport à l'orientation radiale (rad2). Le même phénomène apparaissait également avec la première configuration de sources. Ce résultat n'est pas surprenant car ajouter une composante radiale à un dipôle superficiel ne modifie que très peu le champ magnétique. En revanche, pour une source radiale périphérique (rad2), l'ajout d'une composante tangentielle modifie beaucoup le champ magnétique recueilli, ce qui augmente considérablement le critère à minimiser. Le fait que le moment dipolaire soit légèrement surestimé pour les sources périphériques tangentielles en MEG seule va dans le sens de cette explication.

### 4.8.1.3 Choix d'un critère pour le couplage

Comme dans la première simulation, on observe une amélioration de l'erreur de localisation, de l'erreur d'orientation et de la taille de l'ellipsoïde avec le critère2 par rapport au critère1. En revanche, critère2 introduit un biais de localisation. Il faut donc chercher le compromis le plus favorable entre le critère1 qui, en présence de bruit, fournit des solutions dispersées autour de la bonne solution et le critère 2 qui fournit des solutions moins dispersées et en moyenne plus proches de la bonne solution, mais dont le nuage de dispersion n'est pas centré sur la bonne solution. De plus, l'erreur sur l'amplitude du moment dipolaire est légèrement détériorée par critère2, en particulier pour les dipôles quasi-radiaux, bien que ce soit moins évident que dans la première simulation (figures 4.8 et 4.9). Cette détérioration de l'erreur sur l'amplitude du dipôle s'explique bien, car le critère2 a tendance à multiplier les données par un facteur inférieur à 1 (voir section 4.6.3), l'amplitude du moment dipolaire est donc souvent sous-estimée, et ce d'autant plus que le rapport signal sur bruit est faible. On peut envisager de cumuler la meilleure précision de localisation de critère2 avec la meilleure estimation du moment dipolaire de critère1 en faisant dans un premier temps une localisation avec critère2, en fixant les paramètres de position et d'orientation des sources, et en utilisant avec ces paramètres critère1 pour estimer l'amplitude du moment dipolaire. En revanche, cette technique ne permet pas de compenser le biais de localisation. Comme malgré ce biais, l'erreur de localisation et la taille de l'ellipsoïde sont améliorées, nous avons conservé critère2 pour les simulations suivantes.

### 4.8.1.4 Effet de l'orientation sur la résolution du problème inverse en MEG : apport des modèles réalistes par rapport aux modèles sphériques

Dans cette configuration des sources, nous avons également effectué la résolution des problèmes inverses en MEG seule avec un modèle sphérique. En effet, les simulations effectuées en l'absence de bruit au chapitre 3 montraient qu'un modèle réaliste en MEG permettait de retrouver les sources radiales et l'orientation des sources. Il restait à vérifier si
cet apport des modèles réalistes persistait en présence de bruit, en particulier dans le cas de l'orientation quasi-radiale qui conduit à un signal très faible.

La figure 4.13 compare les erreurs de localisation et d'orientation obtenues avec un modèle réaliste et un modèle sphérique en MEG en présence de bruit. Les dipôles quasi-radiaux ou profonds conduisent en MEG réaliste comme en MEG sphérique à des erreurs de localisation plus importantes que les dipôles à la fois tangentiels et superficiels. Toutefois les erreurs de localisation avec le modèle réaliste restent limitées à une dizaine de millimètres pour les dipôles radiaux, alors que le rapport signal sur bruit est inférieur à 1 (voir figure 4.11), tandis que le modèle sphérique aboutit à des erreurs de localisation de plus de 50 millimètres. En revanche, les erreurs d'orientation obtenues avec un modèle réaliste restent du même ordre, que le dipôle soit radial ou tangentiel, et restent inférieures à 20° ou 30° selon le bruit ajouté. Les erreurs d'orientations obtenues avec le modèle sphérique sont très importantes pour les dipôles radiaux. Pour les dipôles tangentiels, le modèle sphérique conduit à de meilleures reconstructions de l'orientation que le modèle réaliste, puisque le dipôle est contraint à rester tangentiel dans le modèle sphérique alors que son orientation est complètement libre dans le modèle réaliste. Les erreurs d'orientation plus importantes pour les dipôles tangentiels périphériques avec le modèle réaliste s'expliquent par une surestimation de la composante quasi-radiale du dipôle qui affecte peu le signal comparativement aux composantes quasi-tangentielles. Lorsque seules les données MEG sont disponibles, une régularisation comme celle proposée par Fuchs et al. [1998b] devrait permettre d'éviter la surestimation de la composante quasi-radiale, et améliorer l'erreur d'orientation.
<figure><img src="image_68.png" /><figcaption>FIG. 4.13: Erreurs de localisation et d'orientation en MEG réaliste ou sphérique et en présence de bruit. rad1, tan1a et tan1b correspondent à la position de source la plus profonde pour les orientations quasi-radiale et quasi-tangentielles. rad2, tan2a, tan2b correspondent à la position de source la plus périphérique pour les orientations quasi-radiale et quasi-tangentielles.</figcaption></figure>
## 4.9 Troisième configuration : une source dans chaque hémisphère

Nous avons voulu étudier également le cas, fréquent lors de la modélisation de données auditives, où deux sources dipolaires sont présentes simultanément, une dans chaque hémisphère. Un dipôle a été placé dans chaque hémisphère, de manière symétrique, à différentes profondeurs sur des axes temporaux, et selon les trois orientations propres (voir figure 4.14).

<figure><img src="image_69.png" /><figcaption>FIG. 4.14: Troisième configuration des sources : 10 profondeurs * 3 orientations (rouge : quasi-radiales, jaune et bleu : quasi-tangentielles) sur un axe temporal, avec deux dipôles positionnés symétriquement dans chaque hémisphère.</figcaption></figure>

Les réalisations de bruit ont été multipliées par un facteur différent pour chaque source de façon à avoir pour l'EEG à 29 voies un rapport signal sur bruit de 5. Les rapports signal sur bruit ainsi obtenus sont représentés figure 4.15.

<figure><img src="image_70.png" /><figcaption>FIG. 4.15: Rapports signal sur bruit obtenus pour la troisième configuration de sources.</figcaption></figure>
<figure><img src="image_71.png" /><figcaption>FIG. 4.16: Pour la troisième configuration de sources (voir figure 4.14) et en moyenne pour les deux dipôles, de haut en bas : l'erreur de localisation D(mm), l'erreur d'orientation θ(°), l'erreur d'amplitude du moment dipolaire ΔM(%), la taille de l'ellipsoïde (mm³), et le biais de localisation (mm). Les barres d'erreur correspondent à une erreur standard.</figcaption></figure>
### 4.9.1 Résultats

Les résultats de la simulation avec bruit1 et critère2 dans la troisième configuration de source sont présentés figure 4.16. On constate que le couplage permet en moyenne une meilleure localisation des sources, en particulier lorsque les dipôles sont profonds et donc plus proches l'un de l'autre (dipôles symétriques).

* Lorsque les dipôles sont tangentiels et périphériques, la MEG seule donne une bonne précision de localisation, et le couplage donne des résultats comparables, sans détérioration.

* Lorsque les dipôles sont radiaux ou profonds, le couplage donne une erreur de localisation meilleure que la MEG seule ou l'EEG seul.

De même, pour l'erreur d'orientation, la précision sur l'amplitude du moment dipolaire, la taille de l'ellipsoïde et le biais de localisation le couplage donne au moins d'aussi bons résultats que la meilleure des deux modalités.

## 4.10 Quatrième configuration : deux sources dans le même hémisphère

Yvert et al. [2001a] ont montré que plusieurs sources dans le même hémisphère et à faible distance les unes des autres peuvent être présentes simultanément. Afin de savoir si le couplage pourrait permettre de mieux les dissocier, nous avons étudié le cas de deux dipôles actifs simultanément et situés tous deux dans le même hémisphère, avec l'un à une profondeur fixe sur l'axe temporal, et l'autre à une distance variable du premier, toujours sur le même axe (voir figure 4.17).

Les réalisations de bruit ont été multipliées par un facteur différent pour chaque source de façon à avoir pour l'EEG à 29 voies un rapport signal sur bruit de 5 (voir figure 4.18).

### 4.10.1 Résultats

Les résultats de la simulation avec bruit1 et critère2 dans la quatrième configuration de source sont présentés figure 4.19. Le couplage apporte dans la plupart des cas une précision supérieure ou égale à la meilleure des deux modalités, à l'exception du cas où le second dipôle est tangentiel, quand le premier est soit proche soit radial. Dans ce cas, le couplage donne une moins bonne erreur d'orientation que l'EEG et suretime davantage le moment dipolaire. Le biais de localisation n'est généralement pas amélioré par le couplage. En revanche, l'erreur de localisation est presque toujours améliorée.
<figure><img src="image_72.png" /><figcaption>FIG. 4.17: Quatrième configuration des sources : 9 profondeurs * 2 orientations pour le premier dipôle (jaune : quasi-radiales, et noir : quasi-tangentielles) et 1 profondeur * 2 orientations pour le second dipôle (rouge : quasi-radiale, et bleu : quasi-tangentielle).</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>FIG. 4.18: Rapports signal sur bruit obtenus pour la quatrième configuration de sources. De gauche à droite : orientation radiale du premier dipôle et tangentielle du second, orientation tangentielle des deux dipôles, orientation radiale des deux dipôles, et orientation tangentielle du premier dipôle et radiale du second. N.B. : Il y a deux cas de dipôles distants de 5 mm.</figcaption></figure>
<figure><img src="image_74.png" /><figcaption>FIG. 4.19: Pour la quatrième configuration de sources (voir figure 4.17) et en moyenne pour les deux dipôles, de haut en bas : l'erreur de localisation D(mm), l'erreur d'orientation θ(°), l'erreur d'amplitude du moment dipolaire ΔM(%), la taille de l'ellipsoïde (mm³), et le biais de localisation (mm). Les barres d'erreur correspondent à une erreur standard. N.B. : Il y a deux cas de dipôles distants de 5 mm.</figcaption></figure>
## 4.11 Conclusions

Les résultats de cette étude suggèrent que :

* Le couplage apporte toujours une précision supérieure ou égale aux modalités séparées, ce qui confirme les résultats déjà publiés dans le cas de dipôles tangentiels dans une sphère par Stok et al. [1987]; Mosher et al. [1993]; Fuchs et al. [1998b] et dans le cas des sources distribuées par Baillet et al. [1999].

* Les données EEG permettent de retrouver les composantes quasi-radiales des sources dans le couplage, et ce même en ajoutant aux données simulées un bruit réaliste corrélé spatialement.

* L'apport du couplage existe même lorsque peu d'électrodes sont utilisées, bien qu'un nombre plus important d'électrodes permette toujours d'améliorer la précision.

* Lorsque plusieurs sources sont proches, l'apport du couplage devient moins évident même s'il n'est pas nul en termes de localisation.

* Le choix du critère2 pondérant les données en fonction du rapport signal sur bruit pour chaque voie permet une amélioration de 1 à 10 mm de la précision de localisation en MEG réaliste pour les orientations quasi-radiales des dipôles, et une amélioration de 2 à 15° de l'erreur d'orientation pour toutes les orientations. Il conduit à une sous-estimation importante du moment dipolaire qui devrait pouvoir être compensée en utilisant le critère1 pour estimer uniquement l'amplitude du moment dipolaire à la fin de la minimisation.

* Même en présence de bruit, l'utilisation d'un modèle réaliste en MEG permet de retrouver la position des sources quasi-radiales avec une erreur de localisation de l'ordre du centimètre (contre environ 5 cm pour le modèle sphérique), et l'orientation des sources peut être retrouvée avec une erreur inférieure à 20° (contre environ 45° pour le modèle sphérique dans le cas des dipôles quasi-radiaux).
# Chapitre 5

## Conclusions et perspectives
## 5.1 Conclusions

L'objectif de ce travail de thèse était d'étudier les méthodes de localisation des générateurs de l'activité électrique cérébrale en MEG et en EEG dans le but d'améliorer leur précision. Plusieurs approches ont été étudiées dans la littérature dans ce but. Il en ressort que les modèles réalistes avec l'approche problème isolé et des maillages localement affinés permettent une amélioration de la précision de la localisation. Dans cette thèse, trois problématiques ont été abordées : l'évaluation de l'apport d'une interpolation linéaire du potentiel dans la méthode des intégrales de surface en EEG, l'apport des modèles réalistes en MEG, et le couplage des données MEG et EEG.

* L'interpolation linéaire du potentiel dans la méthode des intégrales de surface a été introduite par De Munck [1992] pour l'EEG et Ferguson et al. [1994] pour la MEG. Elle était présentée comme une amélioration par rapport à l'approche centre de masse utilisée jusqu'alors. Bien que plusieurs auteurs l'utilisent [Homma et al., 1994; Schlitt et al., 1995; Zanow, 1997], nous n'avions pu trouver aucune évaluation systématique de l'amélioration apportée par cette méthode. L'étude effectuée dans cette thèse montre que cette approche entraîne une détérioration de la précision pour un même maillage par rapport à l'approche centre de masse, ce qui va dans le même sens que les résultats publiés depuis par Mosher et al. [1999]. En revanche, en choisissant un maillage plus fin conduisant à un temps de calcul identique, on améliore légèrement la précision (de moins d'un millimètre en l'absence de bruit). Plusieurs pistes peuvent être envisagées pour améliorer davantage la précision en modifiant le type de méthode d'intégrales de surfaces : l'augmentation du degré d'interpolation de la surface (triangles courbes) ou l'utilisation d'une formulation variationnelle (méthode de Galerkin) de la méthode des intégrales de surface. Cette approche semble donner de bons résultats [Mosher et al., 1999] mais n'a pas à notre connaissance été évaluée en problème inverse.

* Nous n'avons pu trouver aucune étude systématique du biais de localisation entre le modèle sphérique et les modèles réalistes à trois couches en MEG. De nombreux auteurs jugeaient qu'un modèle réaliste à une couche était suffisant en MEG, sans qu'une étude comparative n'ait pu être trouvée. La plupart des études considéraient des dipôles tangentiels et périphériques, et le comportement des modèles réalistes pour les dipôles d'orientation quelconque était peu décrit. Ces aspects ont été étudiés dans cette thèse, et il en ressort les points suivants. Le biais de localisation avec un modèle sphérique croît en fonction de la position selon un axe vertical de la source (de 2.5 mm en haut de la tête à 12 mm en bas de la tête) et décroît en fonction de la tangentialité (de 10 mm quand le dipôle est radial à 5 mm qui il est tangentiel). De plus, pour un même temps de calcul, un modèle réaliste à trois couches permet une amélioration de 2 mm de
l'erreur de localisation par rapport à un modèle réaliste à une couche. L'orientation des dipôles peut être retrouvée avec un modèle réaliste, et les dipôles radiaux peuvent engendrer des champs magnétiques non négligeables ayant des topographies qui peuvent sembler avoir été engendrées par plusieurs dipôles. Les modèles réalistes, contrairement aux modèles sphériques, conduisent à des erreurs de localisation et d'orientation peu différentes que le dipôle soit radial ou tangentiel quand les données ne sont pas bruitées. Ces résultats ont été obtenus avec des gradiomètres axiaux, ils restent à confirmer dans le cas de gradiomètres planaires. Quand les données sont bruitées, les modèles réalistes retrouvent toujours mieux que les modèles sphériques la position et l'orientation des sources, mais il y a une dégradation pour les dipôles radiaux ou profonds. Fuchs et al. [1998b] proposent d'utiliser une méthode de régularisation pour empêcher que les composantes quasi-radiales des sources soient surestimées. Toutefois aucune comparaison systématique de la précision de la localisation avec ou sans cette régularisation ne semble avoir été effectuée. Wolters et al. [1999] montrent toutefois l'intérêt d'une technique de régularisation similaire pour mieux retrouver le nombre de dipôles pouvant expliquer les données en EEG, et limiter l'amplitude des configurations de sources ne produisant pas de signal. Intuitivement, dans le cas du couplage, cette approche risque de défavoriser la localisation des sources effectivement radiales, tandis que l'ajout de données de type EEG permet d'apporter des informations supplémentaires sur les composantes radiales et les sources profondes sans les pénaliser si elles sont importantes.

* Le couplage de données MEG et EEG était une perspective reconnue pour améliorer la localisation des générateurs cérébraux. Toutefois, dans les études que nous avons pu trouver sur le sujet, soit les deux modalités étaient utilisées séparément et non simultanément, soit les données utilisées avaient été obtenues avec un modèle sphérique. Dans cette thèse, deux méthodes de couplage sont proposées. Elles permettent d'apporter une précision supérieure ou égale à la meilleure des deux modalités séparées, et les composantes quasi-radiales des sources sont mieux retrouvées même en présence de bruit. Cette amélioration est présente même avec un petit nombre d'électrodes dans le couplage (29 électrodes). Toutefois un plus grand nombre d'électrodes permet d'améliorer encore la précision.

## 5.2 Perspectives

Cette thèse propose donc plusieurs techniques d'amélioration de la localisation des sources de l'activité électrique cérébrale. Les études présentées dans ce rapport confirment la difficulté
de la résolution de ce problème inverse, car il existe de nombreuses sources d'erreur et une
bonne précision de localisation ne peut être obtenue que si l'on agit sur chacune de ces causes.

L'une des sources d'erreur qui a été étudiée dans cette thèse est l'erreur sur la définition du
modèle géométrique. Puisque le modèle sphérique entraîne un biais de localisation de
plusieurs millimètres par rapport à un modèle réaliste à trois couches, puisque le modèle à une
couche fait 2 mm d'erreur supplémentaire par rapport au modèle à trois couches, il convient
de décrire avec la plus grande exactitude possible les différents milieux conducteurs qui
constituent la tête. À cet égard, l'idéal serait sans doute d'obtenir une IRM de conductivité de
la tête. Fuchs et al. [1998b] suggèrent une alternative qui consiste à retrouver la conductivité à
partir de données simultanées MEG et EEG, mais la mise en œuvre de cette technique en
pratique est problématique. Si l'ajout de nouvelles inconnues rend le problème encore plus
difficile à résoudre, il persiste aujourd'hui une incertitude sur les valeurs de ces conductivités
qui doivent de plus être différentes selon chaque sujet.

L'imprécision sur la localisation engendrée par l'erreur de modélisation des générateurs cérébraux n'a pas été étudiée dans cette thèse. Toutefois, le modèle dipolaire n'est qu'une approximation du modèle multipolaire. Certains auteurs proposent de revenir à ce modèle [Nolte et Curio, 1997; 2000]. Les modèles de sources distribuées offrent également une alternative intéressante au modèle dipolaire : le problème est linéarisé et la question du nombre de sources ne se pose plus. Toutefois le problème est sous-déterminé, ce qui nécessite l'introduction de nouvelles hypothèses comme par exemple la recherche de la solution de norme minimale, qui ne sont pas toujours justifiables d'un point de vue physiologique.

Comme nous avons montré que l'effet de la position du dipôle selon un axe vertical trouvé
en EEG par Yvert et al. [1997] était également présent en MEG, il est vraisemblable que cet
effet soit également présent dans le couplage, mais il reste à évaluer si cet effet est aussi
important. Ce problème devrait également se retrouver dans le cas de sources distribuées. De
même, l'influence du nombre de couches en MEG et du choix du critère à minimiser devraient
se retrouver pour les sources distribuées.

Les erreurs de localisation sont liées au fait que le problème est mal posé : non seulement
il peut y avoir plusieurs solutions comme c'est le cas des sources radiales en MEG avec un
modèle sphérique mais de plus de petites erreurs sur les calculs ou la définition des paramètres
du modèle peuvent engendrer de grandes erreurs sur la localisation. Lorsqu'un problème
inverse se révèle mal posé, il est habituel d'explorer plusieurs axes de recherche : une autre
formulation du problème permettrait-elle de mieux poser le problème ? Quels sont les
paramètres à estimer ? Les mesures sont-elles suffisantes pour identifier les paramètres que
l'on cherche ? Si le problème reste mal posé, on se ramène à l'optimisation d'une fonction de
coût : laquelle choisir ? Peut-on introduire une technique de régularisation ? Le couplage des
données MEG et EEG peut sans doute être considéré comme une façon de mieux poser le
problème ou d'ajouter des mesures. La prise en compte de contraintes anatomiques (sources situées dans un ruban cortical) ou fonctionnelles (IRMf) pourrait également permettre de mieux poser le problème.

Dans sa thèse, Schwartz [1998] montre que l'erreur de recalage entre l'IRM et les positions des voies MEG ou EEG peut atteindre 3 mm même lorsque l'on procède de façon bien méthodique pour numériser ces positions. De telles erreurs de recalage se répercutent inévitablement sur la précision de la localisation. En revanche, il peut être utile de bruiter les positions individuelles des électrodes et la position de l'ensemble des voies MEG afin de mieux évaluer l'intervalle de confiance des solutions trouvées.

Enfin, une autre source d'erreur est le bruit dans les données. Nous avons vu au chapitre 4 que ce bruit corrélé spatialement peut entraîner une variation de plusieurs millimètres sur la localisation. Il faut donc agir en amont en limitant au maximum les sources de bruit au moment de l'enregistrement, puis ensuite en tenant compte de ce que l'on peut connaître du bruit dans la résolution de notre problème.

Les évaluations réalisées dans cette thèse ont été effectuées par simulation. Cette technique d'évaluation permet une meilleure compréhension de l'influence intrinsèque des différents paramètres de notre modèle : géométrie, conductivité, paramètres spécifiques à la méthode numérique de résolution, etc. Toutefois les simulations ne permettent pas de prédire l'erreur qui sera réalisée en pratique car elles ne peuvent pas complètement rendre compte de toutes les sources de bruit, y compris physiologique, des incertitudes sur les conductivités, de l'erreur dûe au modèle, ou des erreurs de recalage. L'utilisation de fantômes réalistes de la tête [Barth et al., 1986; Menninghaus et al., 1994; Leahy et al., 1998; Baillet et al., 2001] permet de mieux rendre compte de l'influence du bruit ambiant, de l'erreur sur le modèle et des erreurs de recalage. Lorsque c'est possible, les enregistrements de dipôles implantés chez des patients épileptiques permettent de rendre compte de la quasi-totalité des sources d'erreur, à l'exception de l'erreur sur le modèle de source. Toutefois ils ne permettent pas d'analyse systématique des erreurs comme cela est possible avec des simulations numériques. Cependant, de tels enregistrements sont précieux, et leur exploitation devrait permettre de mieux prévoir l'erreur commise par les différentes techniques de localisation des sources de l'activité électromagnétique cérébrale.
<figure><img src="image_75.png" /></figure>
# Annexe A

## Le potentiel et le champ magnétique dans
le modèle sphérique

### 1 Le potentiel électrique [Ary et al., 1981]

Soit un milieu conducteur formé de trois sphères emboîtées, de rayons $R_1, R_2$, et $R_3$, et de conductivités respectives $\sigma_1 = \sigma_{scalp} = \sigma, \sigma_2 = \sigma_{os},$ et $\sigma_3 = \sigma_{cerveau} = \sigma$, et un dipôle d'excentricité $e = z/R_1$ sur l'axe ($Oz$), et dont le moment dipolaire a pour composante radiale $M_r$ et pour composante tangentielle $M_t$.

<figure><img src="image_76.png" /><figcaption>FIG. A.1: Description du modèle sphérique en EEG</figcaption></figure>

<footnote n="0">⁰Cette annexe est référencée aux sections 1.2.1.3, page 15 et 3.1.1, page 58</footnote>
Ary et al. [1981] ont montré que le potentiel électrique mesuré à l'électrode de coordonnées sphériques $(R_1, \theta, \phi)$ s'écrit :

$$V(\theta, \phi) = \frac{1}{4\pi\sigma R_1^2} \sum_{n=1}^{\infty} \frac{2n+1}{n} e^{n-1} \left[ \frac{\zeta(2n+1)^2}{(n+1)d_n} \right] \left[ nM_r P_n(\cos\theta) + M_t P_n^1(\cos\theta) \cos\phi \right]$$

où :

* $P_n$ et $P_n^1$ désignent respectivement le polynôme de Legendre de degré n et son premier polynôme de legendre associé.

* $\zeta = \frac{\sigma_{os}}{\sigma}$

* $d_n = [(n+1)\zeta + n] \left[ \frac{n}{n+1}\zeta + 1 \right] + (1-\zeta) [(n+1)\zeta + n] (f_1^{2n+1} - f_2^{2n+1}) - n(1-\zeta)^2 \left(\frac{f_1}{f_2}\right)^{2n+1}$

* $f_1 = \frac{R_3}{R_1}$ et $f_2 = \frac{R_2}{R_1}$

## 2 Le champ magnétique [Sarvas, 1987]

De même que pour l'EEG, on peut trouver une solution analytique au problème à symétrie sphérique en MEG. Rappelons l'expression du champ magnétique ((1.2.1.11), page 17).

$$\vec{B}(\vec{r}) = \vec{B}_o(\vec{r}) + \frac{\mu_o}{4\pi} \sum_{j=1}^{nbsurf} (\sigma_j^- - \sigma_j^+) \int_{S_j - \exists \vec{r}} V(\vec{r}') \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} \wedge d\vec{S}'$$

Calculons la composante radiale de ce champ magnétique :

$$\vec{B}(\vec{r}) \cdot \frac{\vec{r}}{|\vec{r}|} = \vec{B}_o(\vec{r}) \cdot \frac{\vec{r}}{|\vec{r}|} + \frac{\mu_o}{4\pi} \sum_{j=1}^{nbsurf} (\sigma_j^- - \sigma_j^+) \int_{S_j - \exists \vec{r}} V(\vec{r}') \left( \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} \wedge d\vec{S}' \right) \cdot \frac{\vec{r}}{|\vec{r}|}$$

or

$$\left( \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} \wedge d\vec{S}' \right) \cdot \frac{\vec{r}}{|\vec{r}|} = \left( - \frac{\vec{r} - \vec{r}'}{|\vec{r} - \vec{r}'|^3} \wedge \frac{\vec{r}'}{|\vec{r}'|} \right) \cdot \frac{\vec{r}}{|\vec{r}|} = \vec{0}$$

donc

$$\vec{B}_r(\vec{r}) = \vec{B}(\vec{r}) \cdot \frac{\vec{r}}{|\vec{r}|} = \vec{B}_o(\vec{r}) \cdot \frac{\vec{r}}{|\vec{r}|} = \frac{\mu_o \vec{Q} \wedge \vec{r}_o}{4\pi |\vec{r} - \vec{r}_o|^3} \cdot \frac{\vec{r}}{|\vec{r}|}$$

Par ailleurs, comme à l'extérieur de la tête, $\sigma = 0$, le rotationnel du champ magnétique est nul à l'extérieur de la tête en vertu de la quatrième équation de Maxwell. Il peut donc être dérivé d'un potentiel scalaire magnétique.
$$ \vec{B} = -\mu_o \vec{\nabla} U $$

La troisième équation de Maxwell indique que la divergence du champ magnétique est nulle, donc U est harmonique, et l'on peut choisir U nulle à l'infini.

$$ \begin{cases} \Delta U = 0 \\ U_\infty = 0 \end{cases} $$

Pour trouver U, on utilise l'équation :

$$ U(\vec{r}) = - \int_r^\infty \frac{\partial U}{\partial r} \left( \vec{r}' \cdot \frac{\vec{r}'}{|\vec{r}'|} \right) dr' = \frac{1}{\mu_o} \int_r^\infty \vec{B}_r(\vec{r}') dr' $$

pour obtenir :

$$ \begin{align*} \vec{B}(\vec{r}) &= \vec{\nabla} \left( \frac{\mu_o}{4\pi} \int_r^\infty \frac{\vec{Q} \wedge \vec{r}_o}{|\vec{r} - \vec{r}_o|^3} \cdot \frac{\vec{r}}{|\vec{r}|} \cdot dr' \right) \\ &= \vec{\nabla} \left( \frac{\mu_o}{4\pi} \frac{(\vec{Q} \wedge \vec{r}_o) \cdot \vec{r}'}{F} \right) \\ &= \frac{\mu_o}{4\pi} \frac{F \vec{Q} \wedge \vec{r}_o - \vec{\nabla} F \cdot [(\vec{Q} \wedge \vec{r}_o) \cdot \vec{r}]}{F^2} \end{align*} $$

où

$$ F = |\vec{r} - \vec{r}_o| (|\vec{r}| |\vec{r} - \vec{r}_o| + \vec{r}(\vec{r} - \vec{r}_o)) $$

et

$$ \vec{\nabla} F = \left( \frac{|\vec{r} - \vec{r}_o|^2}{|\vec{r}|} + \frac{(\vec{r} - \vec{r}_o) \cdot \vec{r}}{|\vec{r} - \vec{r}_o|} + 2|\vec{r} - \vec{r}_o| + 2|\vec{r}| \right) \cdot \vec{r} - \left( |\vec{r} - \vec{r}_o| + 2|\vec{r}| + \frac{(\vec{r} - \vec{r}_o) \cdot \vec{r}'}{|\vec{r} - \vec{r}_o|} \right) \cdot \vec{r}_o $$
<figure><img src="image_77.png" /></figure>
# Annexe B

## Équations intégrales du potentiel et du champ magnétique

### 1 Équation intégrale du potentiel

#### 1.1 Résolution de $\Delta u = f$ avec $u \in D$

La méthode des éléments frontières s'appuie sur la connaissance de la résolvante fondamentale du Laplacien en trois dimensions. En effet, au sens des distributions, on vérifie aisément, en prenant l'expression du Laplacien en coordonnées sphériques que :

$$ \Delta \frac{-1}{4\pi|\vec{r} - \vec{r}_o|} = \delta(\vec{r} - \vec{r}_o) $$

Pour résoudre les équations du type $\Delta u = f$ avec $u \in D$, on utilise alors $\langle \Delta u, \phi \rangle = \langle u, \Delta \phi \rangle$, et on l'applique à la résolvante du Laplacien. On obtient ainsi :

$$ u(\vec{r}_o) = \langle u, \delta(\vec{r} - \vec{r}_o) \rangle = \langle \Delta u, \frac{-1}{4\pi|\vec{r} - \vec{r}_o|} \rangle = \langle f, \frac{-1}{4\pi|\vec{r} - \vec{r}_o|} \rangle $$

L'application de cette formule permet le calcul de $u$ en tout point sous forme d'équation intégrale.

#### 1.2 Application à l'équation de Poisson [Geselowitz, 1967]

Dans notre cas, nous cherchons à résoudre le problème (1.2.1.4), page 14 :

$$ \left\{ \begin{array}{l} \nabla(\sigma\vec{\nabla}V) = \nabla(J_p) \\ [V] = 0 \\ [\sigma\vec{\nabla}V \cdot \vec{n}] = 0 \end{array} \right. $$

<footnote n="0">⁰Cette annexe est référencée à la section 1.2.3.1, page 16</footnote>
dans le cas où l'espace G est constitué de N compartiments $G_j$, $j = 1..N$ dans lesquels $\sigma$ est constant. Ces compartiments sont séparés par les N surfaces fermées $S_j$, $j = 1..N$. Pour distinguer la conductivité de part et d'autre de la surface $S_j$, on la note $\sigma_j^+$ à l'extérieur de $S_j$ et $\sigma_j^-$ à l'intérieur.

Nous ne pouvons pas appliquer ce que nous venons de voir à la section précédente, puisque dans notre cas, $u = \sigma V$ présente des discontinuités. L'équation $< \Delta u, \phi > = < u, \Delta \phi >$ est alors à remplacer par la deuxième formule de Green : $< \Delta u, \phi > - < u, \Delta \phi > = [\nabla u \cdot \phi - \nabla \phi \cdot u]$. On obtient alors :

$$< \Delta(\sigma V), \frac{1}{|\vec{r} - \vec{r}_o|} > - < \sigma V, \Delta(\frac{1}{|\vec{r} - \vec{r}_o|}) > = [\nabla(\sigma V) \cdot \frac{1}{|\vec{r} - \vec{r}_o|} - \nabla \frac{1}{|\vec{r} - \vec{r}_o|} \cdot (\sigma V)]$$

Ainsi, en tout point $\vec{r}_o$ non situé sur une surface de discontinuité :

$$< \nabla(\vec{J}_p), \frac{1}{|\vec{r} - \vec{r}_o|} > + 4\pi\sigma(\vec{r}_o)V(\vec{r}_o) = [\nabla(\sigma V) \cdot \frac{1}{|\vec{r} - \vec{r}_o|} - \nabla \frac{1}{|\vec{r} - \vec{r}_o|} \cdot (\sigma V)]$$

Ce qui s'écrit sous forme intégrale :

$$\int_G \nabla(\vec{J}_p) \cdot \frac{1}{|\vec{r} - \vec{r}_o|} dv + 4\pi\sigma(\vec{r}_o)V(\vec{r}_o) = \sum_{j=1}^{N} \int_{S_j} \frac{1}{|\vec{r} - \vec{r}_o|} \cdot \vec{\nabla}(\sigma V) \cdot d\vec{S}_j - \int_{S_j} [(\sigma V) \cdot \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}_o|}] \cdot d\vec{S}_j$$

La continuité du potentiel V, de $1/|\vec{r} - \vec{r}_o|$ et de son gradient, ainsi que le fait que $\sigma$ est constant de part et d'autre des surfaces de discontinuités donne :

$$\int_G \nabla(\vec{J}_p) \cdot \frac{1}{|\vec{r} - \vec{r}_o|} dv + 4\pi\sigma(\vec{r}_o)V(\vec{r}_o) = \sum_{j=1}^{N} \int_{S_j} \frac{1}{|\vec{r} - \vec{r}_o|} \cdot [\sigma \vec{\nabla}(V) \cdot \vec{n}] \cdot dS_j - \int_{S_j} [\sigma] V \cdot \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}_o|} \cdot dS_j$$

Or, les conditions aux limites de notre problème donnent $[\sigma \vec{\nabla}(V) \cdot \vec{n}] = 0$. Donc :

$$\sigma(\vec{r}_o)V(\vec{r}_o) = V_o(\vec{r}_o) - \frac{1}{4\pi} \sum_{j=1}^{N} (\sigma_j^- - \sigma_j^+) \int_{S_j} V \cdot \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}_o|} \cdot d\vec{S}_j \quad (B.1)$$

avec :

$$V_o(\vec{r}_o) = \frac{1}{4\pi} \int_G \vec{J}_p \cdot \vec{\nabla} \left( \frac{1}{|\vec{r} - \vec{r}_o|} \right) dv$$

On note au passage dans cette équation que $V_o$ correspond au potentiel qui serait engendré par les mêmes sources dans un milieu homogène de conductivité unitaire ($\sigma = 1$).

En pratique, pour l'EEG, on mesure $V$ à la surface du scalp, donc sur une surface de discontinuité. Cette formule n'est donc pas applicable directement. Cependant, en tenant compte de la continuité du potentiel $V$, on le calcule par passage à la limite en intégrant l'équation (B.1) sur une petite boule $B_\epsilon$, de rayon $\epsilon$, et centrée au point $\vec{r}_o$ où l'on veut connaître $V$ et appartenant à la surface $S_i$ :

$$\int_{B_\epsilon} \sigma(\vec{r})V(\vec{r})dv = \int_{B_\epsilon} \left( V_o(\vec{r}) - \frac{1}{4\pi} \sum_{j=1}^{N} (\sigma_j^- - \sigma_j^+) \int_{S_j} V \cdot \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}_o|} \cdot d\vec{S}_j \right) dv$$
Or, quand $\epsilon$ tend vers 0, on a les équivalences :

$$ \int_{B_{\epsilon}} \sigma(\vec{r})V(\vec{r})dv \sim \frac{4}{3}\pi\epsilon^3 V(\vec{r}_o) \left( \frac{\Omega_{i,r_o}^{+}}{4\pi}\sigma_i^{+} + \frac{\Omega_{i,r_o}^{-}}{4\pi}\sigma_i^{-} \right) $$

et

$$ \int_{B_{\epsilon}} \left( V_o(\vec{r}) - \frac{1}{4\pi} \sum_{j=1}^{N} (\sigma_j^{-} - \sigma_j^{+}) \int_{S_j} V \cdot \vec{\nabla} \frac{1}{|\vec{r}' - \vec{r}|} \cdot d\vec{S}_j' \right) dv \sim \frac{4}{3}\pi\epsilon^3 \left( V_o(\vec{r}_o) - \frac{1}{4\pi} \sum_{j=1}^{N} (\sigma_j^{-} - \sigma_j^{+}) \int_{S_j \setminus \exists \vec{r}_o} V \cdot \vec{\nabla} \frac{1}{|\vec{r}' - \vec{r}_o|} \cdot d\vec{S}_j' \right) $$

où $\Omega_{i,r_o}^+$ et $\Omega_{i,r_o}^-$ sont les angles solides dont le point $\vec{r}_o$ voit la surface $S_i$ sur laquelle il se trouve. L'intégrale de surface sur $S_j$ n'est pas définie quand $\vec{r}$ se trouve sur $S_j$, c'est pourquoi il faut prendre sa valeur principale, notée $\int_{S_j \setminus \exists \vec{r}_o}$, et qui correspond à la limite de cette intégrale calculée sur la surface $S_j$ privée d'un petit voisinage de $\vec{r}_o$ lorsque la taille de ce voisinage tend vers zéro.

On aboutit ainsi à l'équation intégrale de Fredholm :

$$ (\Omega_{i,r_o}^{+} \sigma_i^{+} + \Omega_{i,r_o}^{-} \sigma_i^{-}) V(\vec{r}_o) = 4\pi V_o(\vec{r}_o) - \sum_{j=1}^{N} (\sigma_j^{-} - \sigma_j^{+}) \int_{S_j \setminus \exists \vec{r}_o} V \cdot \vec{\nabla} \frac{1}{|\vec{r}' - \vec{r}_o|} \cdot d\vec{S}_j' $$

## 2 Équation intégrale du champ magnétique [Geselowitz, 1970]

On part de l'équation (1.2.1.6), page 14 du champ magnétique :

$$ \vec{B}(\vec{r}) = \vec{B}_o(\vec{r}) - \frac{\mu_o}{4\pi} \int_G \sigma \vec{\nabla} V \wedge \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} dv' $$

On utilise l'identité

$$ \vec{\nabla} \wedge (V \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|}) = \vec{\nabla} V \wedge \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} $$

pour obtenir :

$$ \vec{B}(\vec{r}) = \vec{B}_o(\vec{r}) - \frac{\mu_o}{4\pi} \int_G \sigma \vec{\nabla} \wedge (V \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|}) dv' $$

L'espace V se décompose en N compartiments $G_j$ imbriqués les uns dans les autres, tel que $G_1$ est délimité à l'extérieur par la surface $S_1$, $G_2$ est délimité à l'extérieur par la surface $S_2$ et à l'intérieur par $S_1$, et ainsi de suite jusqu'à $S_N$. A l'intérieur de chaque compartiment $G_i$, la conductivité est constante. On la note $\sigma_i^-$, c'est la conductivité à l'intérieur de $S_i$. La
conductivité à l'extérieur de $S_i$ est notée $\sigma_i^+$. Elle est égale à la $\sigma_{i+1}^-$, sauf si $i = N$, où elle vaut 0.

On a alors :

$$ \vec{B}(\vec{r}) = \vec{B}_o(\vec{r}) - \frac{\mu_o}{4\pi} \sum_{i=1}^{N} \sigma_i^- \int_{G_i} \vec{\nabla} \wedge (V\vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|}) dv' $$

Or, si on note $\overline{G_i} = G_1 U \dots U G_i$ le théorème de Gauss donne :

$$ \int_{\overline{G_i}} \vec{\nabla} \wedge (V\vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|}) dv' = - \oint_{S_i} (V\vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|}) \wedge d\vec{S}_i' $$

On utilise alors $\overline{G_1} = G_1$, et $G_i = \overline{G_i} - \overline{G_{i-1}}$ pour $i = 2..N$ :

$$ \vec{B}(\vec{r}) = \vec{B}_o(\vec{r}) + \frac{\mu_o}{4\pi} \sum_{i=1}^{N} \sigma_i^- \oint_{S_i} (V\vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|}) \wedge d\vec{S}_i' \\ - \frac{\mu_o}{4\pi} \sum_{i=2}^{N} \sigma_i^- \oint_{S_{i-1}} (V\vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|}) \wedge d\vec{S}_{i-1}' $$

Un simple changement d'indice en utilisant $\sigma_i^+ = \sigma_{i+1}^-$ et $\sigma_N^+ = 0$ donne alors :

$$ \vec{B}(\vec{r}) = \vec{B}_o(\vec{r}) + \frac{\mu_o}{4\pi} \sum_{i=1}^{N} (\sigma_i^- - \sigma_i^+) \oint_{S_i} (V\vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|}) \wedge d\vec{S}_i' $$
# Annexe C

## Calcul des matrices dans la méthode BEM

### 1 Calcul des matrices BEM pour le potentiel électrique

L'équation linéaire à résoudre est l'équation de Fredholm (1.2.1.9) $(\vec{r} \in S_i)$ :

$$ (\sigma_i^+ \Omega_{i,r}^+ + \sigma_i^- \Omega_{i,r}^-) V(\vec{r}) = 4\pi V_o(\vec{r}) - \sum_{j=1}^{nbsurf} (\sigma_j^- - \sigma_j^+) \int_{S_j - \exists \vec{r}} V(\vec{r}') \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} d\vec{S} $$

Elle est de la forme $L(u) = b$ avec :

$$ \left\{ \begin{aligned} u &= V \\ L(u) &= \frac{\sigma_i^+ \Omega_{i,r}^+ + \sigma_i^- \Omega_{i,r}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} u + \frac{1}{2\pi} \sum_{j=1}^{nbsurf} \frac{\sigma_j^- - \sigma_j^+}{\sigma_i^+ + \sigma_i^-} \int_{S_j - \exists \vec{r}} u(\vec{r}') \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} d\vec{S} \\ b &= \frac{4\pi}{2\pi(\sigma_i^+ + \sigma_i^-)} V_o \end{aligned} \right. $$

L'équation matricielle sera donc de la forme :

$$ AV = V_o $$

où

$$ A = (L(\phi_k)(\vec{x}_j))_{1 \le j, k \le n}, \quad V = (V(\vec{x}_k))_{1 \le k \le n}, \quad V_o = (b(\vec{x}_k))_{1 \le k \le n} $$

#### 1.1 L'approche centre de masse [Hämäläinen, 1987]

Dans cette approche, les nœuds $\vec{r}_k$ sont les centres de gravité des N triangles, notés $\Delta_k$.

$$ \phi_k(\vec{r}) = \begin{cases} 1 & \text{si } \vec{r} \in \Delta_k, \\ 0 & \text{sinon.} \end{cases} $$

<footnote n="0">⁰Cette annexe est référencée aux sections 1.2.3.2, page 19, 1.3.2, page 35, 2.1.2, page 43, et à l'annexe F, page 153</footnote>
On a donc $A_{n,p} = L(\phi_p)(\vec{r}_n)$, soit, si $\vec{r}_n \in S_i$ et $\vec{r}_p \in S_j$ :

$$A_{n,p} = \frac{\sigma_i^+ \Omega_{i,r_n}^+ + \sigma_i^- \Omega_{i,r_n}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} \phi_p(\vec{r}_n) + \frac{1}{2\pi} \sum_{k=1}^{nbsurf} \frac{\sigma_k^- - \sigma_k^+}{\sigma_i^+ + \sigma_i^-} \int_{S_k - \exists \vec{r}_n} \phi_p(\vec{r}) \vec{\nabla} \frac{1}{|\vec{r}_n - \vec{r}|} d\vec{S}$$

Les angles solides dont les centres de gravité des triangles voient l'intérieur et l'extérieur du maillage sur lequel ils se trouvent valent $2\pi$ puisque la surface est régulière en ces points. De plus, on a $\phi_p(\vec{r}_n) = \delta_{p,n}$, et l'intégrale de surface s'annule sur tous les autres triangles que $\Delta_p$. Soit $j$ le numéro de la surface sur laquelle $\vec{r}_p$ se trouve :

$$A_{n,p} = \delta_{p,n} + \frac{1}{2\pi} \frac{\sigma_j^- - \sigma_j^+}{\sigma_i^+ + \sigma_i^-} \int_{\Delta_p - \exists \vec{r}_n} \vec{\nabla} \frac{1}{|\vec{r}_n - \vec{r}|} d\vec{S}$$

On obtient alors l'expression trouvée par Hämäläinen [1987] :

$$A_{n,p} = \delta_{p,n} - \frac{1}{2\pi} \frac{\sigma_j^- - \sigma_j^+}{\sigma_i^+ + \sigma_i^-} \Omega_{p,n}$$

où $\Omega_{p,n}$ est l'angle solide dont le point $\vec{r}_n$ voit le triangle $\Delta_p$.

## 1.2 L'approche interpolation linéaire [De Munck, 1992]

Dans cette approche, les nœuds $\vec{r}_k$ sont les $N$ sommets des triangles du maillage notés $\Delta_{kln}$ si $\vec{r}_k$, $\vec{r}_l$ et $\vec{r}_n$ sont ses sommets.

$$\phi_k(\vec{r}) = \begin{cases} \frac{\det(\vec{r}, \vec{r}_l, \vec{r}_n)}{\det(\vec{r}_k, \vec{r}_l, \vec{r}_n)} & \text{si } \exists l, n, \vec{r} \in \Delta_{kln}, \\ 0 & \text{sinon.} \end{cases}$$

On a donc $A_{m,n} = L(\phi_n)(\vec{r}_m)$, soit, si $\vec{r}_m \in S_i$ et $\vec{r}_n \in \Delta_{kln} \in S_j$ :

$$A_{m,n} = \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} \phi_n(\vec{r}_m) + \frac{1}{2\pi} \sum_{k=1}^{nbsurf} \frac{\sigma_k^- - \sigma_k^+}{\sigma_i^+ + \sigma_i^-} \int_{S_k - \exists \vec{r}_m} \phi_n(\vec{r}) \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} d\vec{S}$$

Or $\phi_n(\vec{r}_m) = \delta_{m,n}$, et l'intégrale de surface est nulle sauf sur les triangles de la surface $S_j$ qui ont pour sommet $\vec{r}_n$ :

$$A_{m,n} = \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} \delta_{m,n} + \frac{1}{2\pi} \frac{\sigma_j^- - \sigma_j^+}{\sigma_i^+ + \sigma_i^-} \sum_{\vec{r}_n \in \Delta_{kln}} \int_{\Delta_{kln} - \exists \vec{r}_m} \phi_n(\vec{r}) \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} d\vec{S}$$

Posons :

$$\beta = \int_{\Delta_{kln} - \exists \vec{r}_m} \phi_n(\vec{r}) \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} d\vec{S}$$
Et faisons le changement d'origine du repère $\vec{y} = \vec{r} - \vec{r}_m$. La formule pour $\phi_n(\vec{y})$ est indépendante de l'origine du repère, on a donc :

$$ \beta = \int_{\Delta_{kln}-\exists \vec{r}_m} \frac{\det(\vec{y}, \vec{y}_k, \vec{y}_l)}{\det(\vec{y}_n, \vec{y}_k, \vec{y}_l)} \vec{\nabla} \frac{1}{|\vec{y}|} \cdot d\vec{S} $$

qui s'écrit également :

$$ \beta = - \int_{\Delta_{kln}-\exists \vec{r}_m} \frac{\det(\vec{y}, \vec{y}_k, \vec{y}_l)}{\det(\vec{y}_n, \vec{y}_k, \vec{y}_l)} \frac{\vec{y}}{|\vec{y}|^3} \cdot d\vec{S} $$

Posons alors $\vec{a} = (\vec{y}_k - \vec{y}_n) \wedge (\vec{y}_l - \vec{y}_n)$. Il s'agit en fait du vecteur normal au triangle dont la norme $|\vec{a}| = A$ vaut deux fois l'aire du triangle. On a alors :

$$ \beta = - \int_{\Delta_{kln}-\exists \vec{r}_m} \frac{\det(\vec{y}, \vec{y}_k, \vec{y}_l)}{\det(\vec{y}_n, \vec{y}_k, \vec{y}_l)} \frac{\vec{y}}{|\vec{y}|^3} \cdot \frac{\vec{a}}{A} dS $$

Or, puisque $\vec{y}_m = \vec{0}$ et $\vec{y} \in \Delta_{kln}$ :

$$ \vec{y} \cdot \frac{\vec{a}}{A} = \det(\vec{y}_m, \Delta_{kln}) = \vec{y}_n \cdot \frac{(\vec{y}_k - \vec{y}_n) \wedge (\vec{y}_l - \vec{y}_n)}{A} = \vec{y}_n \cdot \frac{\vec{y}_k \wedge \vec{y}_l}{A} = \frac{\det(\vec{y}_n, \vec{y}_k, \vec{y}_l)}{A} $$

D'où :

$$ \beta = -\frac{1}{A} \int_{\Delta_{kln}-\exists \vec{r}_m} \frac{\det(\vec{y}, \vec{y}_k, \vec{y}_l)}{|\vec{y}|^3} dS $$

Et :

$$ \beta = -\frac{1}{A} \int_{\Delta_{kln}-\exists \vec{r}_m} \frac{\vec{y} \cdot (\vec{y}_k \wedge \vec{y}_l)}{|\vec{y}|^3} dS $$

Pour simplifier l'écriture des calculs, nous allons changer de notation. Il nous faut calculer, dans un triangle quelconque $\Delta_{123}$, la valeur des :

$$ \Omega_i = \frac{1}{A} \int_{\Delta_{123}} \frac{\vec{y} \cdot \vec{z}_i}{|\vec{y}|^3} dS, \quad i = 1, 2, 3. $$

lorsque $\vec{z}_i = \vec{y}_{i+1} \wedge \vec{y}_{i-1}$ avec, par convention, des indices définis modulo 3 : $\vec{y}_0 = \vec{y}_3$, et $\vec{y}_1 = \vec{y}_4$.

Pour calculer les $\Omega_i$, on introduit le vecteur $\vec{\Omega}$ défini par :

$$ \vec{\Omega} = \int_{\Delta_{123}} \vec{\nabla} \left( \frac{1}{|\vec{y}|} \right) \wedge d\vec{S} \qquad (C.1) $$

On a alors :

$$ \vec{\Omega} = - \oint_{\partial \Delta_{123}} \frac{d\vec{\gamma}}{|\vec{\gamma}|} = - \sum_{i=1}^{3} \int_{0}^{1} \frac{\vec{y}_{i+1} - \vec{y}_{i}}{|\vec{y}_{i} + s(\vec{y}_{i+1} - \vec{y}_{i})|} ds $$
que l'on peut intégrer explicitement :

$$
\vec{\Omega} = \sum_{i=1}^{3} \left\{ \begin{array}{ll} \left[ \frac{\ln((\vec{y}_i + s(\vec{y}_{i+1} - \vec{y}_i))(\vec{y}_{i+1} - \vec{y}_i) + |\vec{y}_i + s(\vec{y}_{i+1} - \vec{y}_i)||\vec{y}_{i+1} - \vec{y}_i|)}{|\vec{y}_{i+1} - \vec{y}_i|} \right]^1_0 & (\vec{y}_{i+1} - \vec{y}_i) \\ \\ \left[ \frac{\ln(|\vec{y}_i + s(\vec{y}_{i+1} - \vec{y}_i)|)}{|\vec{y}_{i+1} - \vec{y}_i|} \right]^1_0 & \text{sinon} \end{array} \right. \\
\text{d'où :}
$$

d'où :

$$
\vec{\Omega} = \sum_{i=1}^{3} \begin{cases} \frac{\vec{y}_{i+1} - \vec{y}_i}{|\vec{y}_{i+1} - \vec{y}_i|} \ln\left(\frac{\vec{y}_{i+1} \cdot (\vec{y}_{i+1} - \vec{y}_i) + |\vec{y}_{i+1}||\vec{y}_{i+1} - \vec{y}_i|}{\vec{y}_i \cdot (\vec{y}_{i+1} - \vec{y}_i) + |\vec{y}_i||\vec{y}_{i+1} - \vec{y}_i|}\right) & \text{si } (\vec{y}_i \cdot (\vec{y}_{i+1} - \vec{y}_i))^2 \neq |\vec{y}_i|^2 |\vec{y}_{i+1} - \vec{y}_i|^2 \\ \\ \frac{\vec{y}_{i+1} - \vec{y}_i}{|\vec{y}_{i+1} - \vec{y}_i|} \ln\left(\frac{|\vec{y}_{i+1}|}{|\vec{y}_i|}\right) & \text{sinon} \end{cases}
$$

Il suffit alors de séparer la somme en deux et de faire un changement d'indice pour obtenir :

$$
\vec{\Omega} = \sum_{i=1}^{3} (\nu_{i-1}^{0} - \nu_{i}^{0}) \vec{y}_{i} \quad (C.2)
$$

où :

$$
\nu_i^0 = \begin{cases} \frac{-1}{|\vec{y}_{i+1} - \vec{y}_i|} \ln\left(\frac{\vec{y}_i \cdot (\vec{y}_{i+1} - \vec{y}_i) + |\vec{y}_i||\vec{y}_{i+1} - \vec{y}_i|}{\vec{y}_{i+1} \cdot (\vec{y}_{i+1} - \vec{y}_i) + |\vec{y}_{i+1}||\vec{y}_{i+1} - \vec{y}_i|}\right) & \text{si } (\vec{y}_i \cdot (\vec{y}_{i+1} - \vec{y}_i))^2 \neq |\vec{y}_i|^2 |\vec{y}_{i+1} - \vec{y}_i|^2 \\ \\ \frac{-1}{|\vec{y}_{i+1} - \vec{y}_i|} \ln\left(\frac{|\vec{y}_i|}{|\vec{y}_{i+1}|}\right) & \text{sinon} \end{cases}
$$

Nous allons à présent montrer le lien entre $\vec{\Omega}$ et les $\Omega_i$ que nous devions calculer. Pour cela, on introduit les coordonnées barycentriques de $\vec{y}$ dans le triangle $\Delta_{123}: \vec{y} = \sum_{j=1}^{3} \lambda_j \vec{y}_j$, $\sum_{j=1}^{3} \lambda_j = 1$. On peut alors exprimer explicitement les coordonnées barycentriques, puisque :

$$
\vec{y} \cdot \vec{z}_i = \sum_{j=1}^{3} \lambda_j \vec{y}_j \cdot \vec{z}_i = \lambda_i \vec{y}_i \cdot \vec{z}_i = \lambda_i \vec{y}_i \cdot (\vec{y}_{i+1} \wedge \vec{y}_{i-1}) = \lambda_i \det(\vec{y}_i, \vec{y}_{i+1}, \vec{y}_{i-1})
$$

d'où :

$$
\lambda_i = \frac{\vec{y} \cdot \vec{z}_i}{\det(\vec{y}_1, \vec{y}_2, \vec{y}_3)}
$$

On injecte alors dans $\vec{\Omega}$ l'expression de $\vec{y}$ en coordonnées barycentriques :

$$
\vec{\Omega} = - \int_{\Delta_{123}} \frac{\vec{y}}{|\vec{y}|^3} \wedge d\vec{S} = - \frac{1}{\det(\vec{y}_1, \vec{y}_2, \vec{y}_3)} \sum_{i=1}^{3} \int_{\Delta_{123}} \vec{y} \cdot \vec{z}_i \frac{\vec{y}_i}{|\vec{y}|^3} \wedge d\vec{S}
$$

d'où :

$$
-\det(\vec{y}_1, \vec{y}_2, \vec{y}_3)\vec{\Omega} = \sum_{i=1}^{3} \int_{\Delta_{123}} \vec{y} \cdot \vec{z}_i \frac{\vec{y}_i}{|\vec{y}|^3} \wedge \frac{\vec{a}}{A} dS = \sum_{i=1}^{3} (\Omega_i \vec{y}_i) \wedge \vec{a}
$$

Ce qui donne un système de trois équations à trois inconnues, mais ne permet pas de trouver
les $\Omega_i$, car ce système est lié : $\det(\vec{y}_1 \wedge \vec{a}, \vec{y}_2 \wedge \vec{a}, \vec{y}_3 \wedge \vec{a}) = 0$. Il nous faut donc une équation
supplémentaire, qui nous est donnée par :

$$
\sum_{i=1}^{3} \Omega_i = \frac{1}{A} \int_{\Delta_{123}} \frac{\vec{y} \cdot (\sum_{i=1}^{3} \vec{z}_i)}{|\vec{y}|^3} dS = \frac{1}{A} \int_{\Delta_{123}} \frac{\vec{y} \cdot \vec{a}}{|\vec{y}|^3} dS = \int_{\Delta_{123}} \frac{\vec{y} \cdot d\vec{S}}{|\vec{y}|^3} = \Omega
$$
où $\Omega$ est l'angle solide dont le point origine $\vec{y}_m$ voit le triangle $\Delta_{123}$. La formule analytique de cet angle solide est connue :

$$ \Omega = 2 \arctan \left( \frac{\det(\vec{y}_1, \vec{y}_2, \vec{y}_3)}{|\vec{y}_1||\vec{y}_2||\vec{y}_3| + |\vec{y}_1|(\vec{y}_2 \cdot \vec{y}_3) + |\vec{y}_2|(\vec{y}_3 \cdot \vec{y}_1) + |\vec{y}_3|(\vec{y}_1 \cdot \vec{y}_2)} \right) $$

Le système à résoudre est donc :

$$ \left\{ \begin{array}{l} -\det(\vec{y}_1, \vec{y}_2, \vec{y}_3)\vec{\Omega} = \sum_{i=1}^{3} (\Omega_i \vec{y}_i) \wedge \vec{a} \\ \sum_{i=1}^{3} \Omega_i = \Omega \end{array} \right. $$

d'où :

$$ \left\{ \begin{array}{l} -\det(\vec{y}_1, \vec{y}_2, \vec{y}_3)\vec{\Omega} \cdot (\vec{y}_{i+1} - \vec{y}_{i-1}) = \sum_{j=1}^{3} \Omega_j (\vec{y}_j \wedge \vec{a}) \cdot (\vec{y}_{i+1} - \vec{y}_{i-1}) \\ \sum_{i=1}^{3} \Omega_i = \Omega \end{array} \right. $$

et :

$$ \left\{ \begin{array}{l} -\det(\vec{y}_1, \vec{y}_2, \vec{y}_3)\vec{\Omega} \cdot (\vec{y}_{i+1} - \vec{y}_{i-1}) = \sum_{j=1}^{3} \Omega_j ((\vec{y}_{i+1} - \vec{y}_{i-1}) \wedge \vec{y}_j) \cdot \vec{a} \\ \sum_{i=1}^{3} \Omega_i = \Omega \end{array} \right. $$

donc :

$$ \left\{ \begin{array}{l} -\det(\vec{y}_1, \vec{y}_2, \vec{y}_3)\vec{\Omega} \cdot (\vec{y}_{i+1} - \vec{y}_{i-1}) = (\Omega_i \vec{z}_{i-1} + \Omega_i \vec{z}_{i+1} - \Omega_{i+1} \vec{z}_i - \Omega_{i-1} \vec{z}_i) \cdot \vec{a} \\ \sum_{i=1}^{3} \Omega_i = \Omega \end{array} \right. $$

On injecte alors la seconde équation dans la première en utilisant $\sum_{i=1}^{3} \vec{z}_i = \vec{a}$ pour trouver :

$$ -\det(\vec{y}_1, \vec{y}_2, \vec{y}_3)\vec{\Omega} \cdot (\vec{y}_{i+1} - \vec{y}_{i-1}) = \Omega_i \vec{a} \cdot \vec{a} - \Omega \vec{z}_i \cdot \vec{a} $$

D'où :

$$ \Omega_i = \frac{1}{A^2} \left( \vec{z}_i \cdot \vec{a} \Omega - \det(\vec{y}_1, \vec{y}_2, \vec{y}_3) (\vec{y}_{i+1} - \vec{y}_{i-1}) \cdot \vec{\Omega} \right) $$

Nous avons alors tous les éléments pour calculer les coefficients de la matrice :

$$ A_{m,n} = \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} \delta_{m,n} - \frac{1}{2\pi} \frac{\sigma_j^- - \sigma_j^+}{\sigma_i^+ + \sigma_i^-} \sum_{\vec{r}_n \in \Delta_{kln}} \Omega_n(\Delta_{kln}, \vec{r}_m) $$

Les termes diagonaux nécessitent l'évaluation des angles solides $\Omega_{i,r_m}^+$ et $\Omega_{i,r_m}^-$, pour lesquels le temps de calcul est assez long. Pour calculer ces termes, on utilisera plutôt une propriété de la matrice, liée au fait que le potentiel est défini à une constante près : le vecteur composé uniquement de 1 est un vecteur propre de la matrice. Montrons le :

$$ \begin{aligned} \sum_{n=1}^{N} A_{m,n} &= \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} + \sum_{n=1}^{N} \frac{1}{2\pi} \sum_{k=1}^{nbsurf} \frac{\sigma_k^- - \sigma_k^+}{\sigma_i^+ + \sigma_i^-} \int_{S_k - \exists \vec{r}_m} \phi_n(\vec{r}) \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} \cdot d\vec{S} \\ &= \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} + \frac{1}{2\pi} \sum_{k=1}^{nbsurf} \frac{\sigma_k^- - \sigma_k^+}{\sigma_i^+ + \sigma_i^-} \int_{S_k - \exists \vec{r}_m} \left( \sum_{n=1}^{N} \phi_n(\vec{r}) \right) \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} \cdot d\vec{S} \end{aligned} $$
Or on vérifie aisément que sur tout triangle du maillage, $\sum_{n=1}^{N} \phi_n(\vec{r}) = 1$, d'où :

$$ 
\begin{align*} 
\sum_{n=1}^{N} A_{m,n} &= \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} + \frac{1}{2\pi} \sum_{k=1}^{nbsurf} \frac{\sigma_k^- - \sigma_k^+}{\sigma_i^+ + \sigma_i^-} \int_{S_k - \exists \vec{r}_m} \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} \cdot d\vec{S} \\
&= \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} - \frac{1}{2\pi} \sum_{k=1}^{nbsurf} \frac{\sigma_k^- - \sigma_k^+}{\sigma_i^+ + \sigma_i^-} \Omega_{k,r_m}^- \\
&= \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^- - (\sigma_i^- - \sigma_i^+) \Omega_{i,r_m}^- - \sum_{k=i+1}^{nbsurf} (\sigma_k^- - \sigma_k^+) 4\pi}{2\pi(\sigma_i^+ + \sigma_i^-)}
\end{align*} 
$$

Or $\sigma_k^- = \sigma_{k-1}^+$, et $\sigma_{nbsurf}^+ = 0$, d'où :

$$ \sum_{n=1}^{N} A_{m,n} = \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^- - (\sigma_i^- - \sigma_i^+) \Omega_{i,r_m}^- - \sigma_i^+ 4\pi}{2\pi(\sigma_i^+ + \sigma_i^-)} $$

$$ \sum_{n=1}^{N} A_{m,n} = \frac{\sigma_i^+(\Omega_{i,r_m}^+ + \Omega_{i,r_m}^- - 4\pi)}{2\pi(\sigma_i^+ + \sigma_i^-)} = 0 $$

Pour calculer les termes diagonaux, il suffit donc d'utiliser l'expression :

$$ A_{m,m} = - \sum_{n=1, n \neq m}^{N} A_{m,n} $$

## 2 Calcul des matrices BEM pour le champ magnétique

Le champ magnétique dépend linéairement du potentiel selon l'équation de Geselowitz (1.2.1.11) :

$$ \vec{B}(\vec{r}) = \vec{B}_o(\vec{r}) + \frac{\mu_o}{4\pi} \sum_{j=1}^{nbsurf} (\sigma_j^- - \sigma_j^+) \int_{S_j - \exists \vec{r}} V(\vec{r}') \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} \wedge d\vec{S}' $$

Pour calculer les composantes du champ magnétique aux M spires de position $\vec{r}_m$ et d'orientation $\vec{n}_m$, on calcule donc préalablement le potentiel électrique aux N nœuds $\vec{r}_n$ du maillage. Comme en tout point des surfaces de discontinuité $S_j$, $V(\vec{r}) = \sum_{i=1}^{N} V(\vec{r}_i)\phi_n(\vec{r})$, notre équation peut s'écrire sous forme matricielle :

$$ B = B_o + MV $$

où :

$$ 
\left\{ 
\begin{aligned}
B &= (\vec{B}(\vec{r}_m) \cdot \vec{n}_m)_{1 \le m \le M} \\
B_o &= (\vec{B}_o(\vec{r}_m) \cdot \vec{n}_m)_{1 \le m \le M} \\
M &= \left( \frac{\mu_o}{4\pi} \sum_{j=1}^{nbsurf} (\sigma_j^- - \sigma_j^+) \int_{S_j - \exists \vec{r}_m} \phi_n(\vec{r}) \left( \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} \wedge d\vec{S} \right) \cdot \vec{n}_m \right)_{1 \le n \le N, 1 \le m \le M} \\
V &= (V(\vec{r}_n))_{1 \le n \le N}
\end{aligned}
\right. 
$$
## 2.1 Approche centre de masse

Dans cette approche, les nœuds $\vec{r}_k$ sont les centres de gravité des N triangles, notés $\Delta_k$.

$$ \phi_k(\vec{r}) = \begin{cases} 1 & \text{si } \vec{r} \in \Delta_k, \\ 0 & \text{sinon.} \end{cases} $$

On a donc, pour $\vec{r}_n \in S_j$ :

$$ M_{m,n} = \frac{\mu_o}{4\pi} (\sigma_j^{-} - \sigma_j^{+}) \int_{\Delta_n - \exists \vec{r}_m} (\vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} \wedge d\vec{S}) \cdot \vec{n}_m $$

Où l'on retrouve le vecteur $\vec{\Omega}$ défini en (C.1) à la section précédente et que l'on peut calculer avec la formule (C.2) :

$$ M_{m,n} = \frac{\mu_o}{4\pi} (\sigma_j^{-} - \sigma_j^{+}) \vec{\Omega}(\Delta_n, \vec{r}_m) \cdot \vec{n}_m $$

## 2.2 Approche interpolation linéaire [Ferguson et al., 1994]

Dans cette approche, les nœuds $\vec{r}_k$ sont les N sommets des triangles du maillage notés $\Delta_{kln}$ si $\vec{r}_k, \vec{r}_l$ et $\vec{r}_n$ sont ses sommets.

$$ \phi_k(\vec{r}) = \begin{cases} \frac{\det(\vec{r}, \vec{r}_l, \vec{r}_n)}{\det(\vec{r}_k, \vec{r}_l, \vec{r}_n)} & \text{si } \exists l, n, \vec{r} \in \Delta_{kln}, \\ 0 & \text{sinon.} \end{cases} $$

On a donc, pour $\vec{r}_n \in S_j$ :

$$ M_{m,n} = \frac{\mu_o}{4\pi} (\sigma_j^{-} - \sigma_j^{+}) \sum_{\vec{r}_n \in \Delta_{kln}} \int_{\Delta_{kln} - \exists \vec{r}_m} \phi_n(\vec{r}) (\vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} \wedge d\vec{S}) \cdot \vec{n}_m $$

On est donc amené à calculer la valeur du vecteur :

$$ \vec{a}_{kln}^{m} = \int_{\Delta_n - \exists \vec{r}_m} \phi_n(\vec{r}) (\vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} \wedge d\vec{S}) $$

On utilise une combinaison du théorème de Stokes :

$$ \int_{\Delta} \vec{\nabla} u \wedge d\vec{S} = - \oint_{\partial \Delta} u d\vec{l} $$

et de

$$ \vec{\nabla}(uv) = u\vec{\nabla}v + v\vec{\nabla}u $$

pour obtenir :

$$ \vec{a}_{kln}^{m} = - \oint_{\partial \Delta_{kln}} \frac{\phi_n(\vec{r})}{|\vec{r}_m - \vec{r}|} d\vec{l} - \int_{\Delta_{kln} - \exists \vec{r}_m} \frac{1}{|\vec{r}_m - \vec{r}|} (\vec{\nabla} \phi_n(\vec{r}) \wedge d\vec{S}) $$
La somme sur tous les triangles du premier terme s'annule car chaque arête des contours des triangles est parcourue une fois dans un sens et une fois dans l'autre. Pour calculer le second terme, on évalue le gradient des fonctions-test pour $\vec{r} \in \Delta_{kln}$ :

$$ \vec{\nabla} \phi_n(\vec{r}) = \vec{\nabla} \frac{\det(\vec{r}, \vec{r}_k, \vec{r}_l)}{\det(\vec{r}_k, \vec{r}_l, \vec{r}_n)} = \vec{\nabla} \frac{\vec{r} \cdot (\vec{r}_k \wedge \vec{r}_l)}{\det(\vec{r}_k, \vec{r}_l, \vec{r}_n)} = \frac{\vec{r}_k \wedge \vec{r}_l}{\det(\vec{r}_k, \vec{r}_l, \vec{r}_n)} $$

Par ailleurs on peut écrire :

$$ d\vec{S} = \frac{(\vec{r}_k - \vec{r}_n) \wedge (\vec{r}_l - \vec{r}_n)}{|(\vec{r}_k - \vec{r}_n) \wedge (\vec{r}_l - \vec{r}_n)|} dS = \frac{(\vec{r}_k - \vec{r}_n) \wedge (\vec{r}_l - \vec{r}_n)}{A} dS $$

où A vaut deux fois l'aire du triangle. On a alors :

$$ \begin{align*} \vec{\nabla} \phi_n(\vec{r}) \wedge d\vec{S} &= \frac{\vec{r}_k \wedge \vec{r}_l}{\det(\vec{r}_k, \vec{r}_l, \vec{r}_n)} \wedge \frac{(\vec{r}_k - \vec{r}_n) \wedge (\vec{r}_l - \vec{r}_n)}{A} dS \\ &= \frac{-(((\vec{r}_k - \vec{r}_n) \wedge (\vec{r}_l - \vec{r}_n)) \cdot \vec{r}_l)\vec{r}_k + (((\vec{r}_k - \vec{r}_n) \wedge (\vec{r}_l - \vec{r}_n)) \cdot \vec{r}_k)\vec{r}_l}{\det(\vec{r}_k, \vec{r}_l, \vec{r}_n)A} dS \\ &= \frac{-\det(\vec{r}_k, \vec{r}_l, \vec{r}_n)\vec{r}_k + \det(\vec{r}_k, \vec{r}_l, \vec{r}_n)\vec{r}_l}{\det(\vec{r}_k, \vec{r}_l, \vec{r}_n)A} dS \\ &= \frac{\vec{r}_l - \vec{r}_k}{A} dS \end{align*} $$

Et finalement :

$$ \vec{a}_{kln}^{m} = - \oint_{\partial \Delta_{kln}} \frac{\phi_n(\vec{r})}{|\vec{r}_m - \vec{r}|} d\vec{l} - \frac{\vec{r}_l - \vec{r}_k}{A} \int_{\Delta_{kln} - \exists \vec{r}_m} \frac{dS}{|\vec{r}_m - \vec{r}|} $$

Il nous faut donc calculer $\int_{\Delta_{kln} - \exists \vec{r}_m} \frac{dS}{|\vec{r}_m - \vec{r}|}$. Pour cela on commence par faire le changement de variable $\vec{y} = \vec{r} - \vec{r}_m$, puis on note A la projection du point $\vec{r}_m$ sur le plan du triangle, et $\vec{c}$ le vecteur ayant pour origine $\vec{r}_m$ et pour extrémité A.

<figure><img src="image_78.png" /></figure>

Pour simplifier les écritures, on peut remplacer les indices $k, l, n$ par 1, 2 et 3 respectivement, en définissant modulo 3 l'indice des points correspondants : $\vec{r}_{i+3} = \vec{r}_i$. On peut alors remarquer que l'intégrale sur le triangle $\Delta_{kln}$ est égale à la somme des intégrales sur les triangles $\Delta_{A,i,i+1}$, pondérée par un signe dépendant de l'orientation du triangle, dont la normale est notée $\vec{n}$.

$$ \int_{\Delta_{kln}} \frac{dS}{|\vec{y}|} = \sum_{i=1}^{3} \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{n})) \int_{\Delta_{A,i,i+1}} \frac{dS}{|\vec{y}|} $$
et on peut écrire :

$$ \int_{\Delta_{A,i,i+1}} \frac{dS}{|\vec{y}|} = \int_{0}^{\theta_{i,i+1}} \int_{0}^{|\vec{c}_{\theta}|} \frac{r dr d\theta}{\sqrt{|\vec{c}|^2 + r^2}} = \int_{0}^{\theta_{i,i+1}} |\vec{y}_{\theta}| d\theta + |\vec{c}| \theta_{i,i+1} $$

Pour évaluer l'angle (strictement positif) $\theta_{i,i+1}$, on introduit l'angle solide sous lequel le triangle $\Delta_{A,i,i+1}$ est vu depuis le point $\vec{r}_m$ :

$$ \begin{aligned} \Omega_{A,i,i+1} &= \int_{\Delta_{A,i,i+1}} \frac{\vec{y} \cdot d\vec{S}}{|\vec{y}|^3} = -|\vec{c}| \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{c})) \int_{0}^{\theta_{i,i+1}} \int_{0}^{|\vec{c}_{\theta}|} \frac{r dr d\theta}{\sqrt{|\vec{c}|^2 + r^2}} \\ &= -|\vec{c}| \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{c})) \int_{0}^{\theta_{i,i+1}} \left[ \frac{1}{\sqrt{|\vec{c}|^2 + r^2}} \right]_{0}^{|\vec{c}_{\theta}|} d\theta \\ &= \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{c})) \theta_{i,i+1} - |\vec{c}| \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{c})) \int_{0}^{\theta_{i,i+1}} \frac{1}{|\vec{y}_{\theta}|} d\theta \end{aligned} $$

En remplaçant $\theta$ dans notre équation, on trouve alors :

$$ \begin{aligned} \int_{\Delta_{A,i,i+1}} \frac{dS}{|\vec{y}|} &= \int_{0}^{\theta_{i,i+1}} \frac{|\vec{y}_{\theta}|^2 - |\vec{c}|^2}{|\vec{y}_{\theta}|} d\theta - |\vec{c}| \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{c})) \Omega_{A,i,i+1} \\ &= \int_{0}^{\theta_{i,i+1}} \frac{|\vec{c}_{\theta}|^2}{|\vec{y}_{\theta}|} d\theta - |\vec{c}| \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{c})) \Omega_{A,i,i+1} \end{aligned} $$

On fait alors le changement de variable $\alpha$ tel que $\vec{c}_{\theta} = \vec{c}_i + \alpha(\vec{c}_{i+1} - \vec{c}_i)$ :

$$ \frac{\cos \theta}{\sin \theta} = \frac{\vec{c}_i \cdot \vec{c}_{\theta}}{|\vec{c}_i \wedge \vec{c}_{\theta}|} $$

d'où, comme $\vec{c}_i \wedge \vec{c}_{\theta} = \alpha \vec{c}_i \wedge \vec{c}_{i+1}$ :

$$ \frac{-1}{\sin^2 \theta} d\theta = \frac{\vec{c}_i \cdot (\vec{c}_{i+1} - \vec{c}_i) |\vec{c}_i \wedge \vec{c}_{\theta}| - \vec{c}_i \cdot \vec{c}_{\theta} |\vec{c}_i \wedge \vec{c}_{i+1}|}{|\vec{c}_i \wedge \vec{c}_{\theta}|^2} d\alpha $$

et comme $|\vec{c}_i \wedge \vec{c}_{\theta}|^2 = |\vec{c}_i|^2 |\vec{c}_{\theta}|^2 \sin^2 \theta$ :

$$ -|\vec{c}_i|^2 |\vec{c}_{\theta}|^2 d\theta = \vec{c}_i \cdot (\alpha(\vec{c}_{i+1} - \vec{c}_i) - \vec{c}_{\theta}) |\vec{c}_i \wedge \vec{c}_{i+1}| d\alpha = -|\vec{c}_i|^2 |\vec{c}_i \wedge \vec{c}_{i+1}| d\alpha $$

soit :

$$ |\vec{c}_{\theta}|^2 d\theta = |\vec{c}_i \wedge \vec{c}_{i+1}| d\alpha $$

On peut alors écrire :

$$ \int_{\Delta_{A,i,i+1}} \frac{dS}{|\vec{y}|} = \int_{0}^{1} \frac{|\vec{c}_i \wedge \vec{c}_{i+1}|}{|\vec{y}_i + \alpha(\vec{y}_{i+1} - \vec{y}_i)|} d\alpha - |\vec{c}| \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{c})) \Omega_{A,i,i+1} $$

d'où :

$$ \int_{\Delta_{A,i,i+1}} \frac{dS}{|\vec{y}|} = |\vec{c}_i \wedge \vec{c}_{i+1}| \nu_i^o - |\vec{c}| \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{c})) \Omega_{A,i,i+1} $$
où l'on retrouve les $\nu_i^o$ définis en (C.2). En regroupant ces résultats, on trouve alors :

$$ 
\begin{align*}
\int_{\Delta_{kln}} \frac{dS}{|\vec{y}|} &= \sum_{i=1}^{3} \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{n})) (|\vec{c}_i \wedge \vec{c}_{i+1}| \nu_i^o - |\vec{c}| \operatorname{sgn}(\det(\vec{c}_i, \vec{c}_{i+1}, \vec{c})) \Omega_{A,i,i+1}) \\
&= \sum_{i=1}^{3} (\vec{c}_i \wedge \vec{c}_{i+1}) \cdot \vec{n} \nu_i^o - \vec{c} \cdot \vec{n} \Omega_{kln}
\end{align*}
$$

et donc :

$$ M_{m,n} = \frac{\mu_o}{4\pi} (\sigma_j^- - \sigma_j^+) \sum_{\vec{r}_n \in \Delta_{kln}} \frac{\vec{r}_k - \vec{r}_l}{A_{\Delta_{kln}}} \left( \sum_{i=1}^{3} (\vec{c}_i \wedge \vec{c}_{i+1}) \cdot \vec{n}_{\Delta_{kln}} \nu_i^o - \vec{c} \cdot \vec{n}_{\Delta_{kln}} \Omega_{kln} \right) $$
# Annexe D

## Influence de la conductivité en EEG et en MEG

Cette annexe montre que l'EEG dépend inversement d'un facteur multiplicatif affectant toutes les conductivités alors que la MEG ne dépend que des rapports entre les conductivités des différents milieux.

### 1 Effet d'un facteur multiplicatif affectant toutes les conductivités en EEG

Le potentiel électrique est donné par l'équation de Fredholm (1.2.1.9 page 16) :

$$ (\sigma_i^+ \Omega_{i,r}^+ + \sigma_i^- \Omega_{i,r}^-) V(\vec{r}) = 4\pi V_o(\vec{r}) - \sum_{j=1}^{nbsurf} (\sigma_j^- - \sigma_j^+) \int_{S_j - \exists \vec{r}} V(\vec{r}') \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} \cdot d\vec{S} $$

où

$$ V_o(\vec{r}) = \frac{1}{4\pi} \int_G \vec{J}_p \cdot \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} dv $$

Le potentiel $V_o$ ne dépend pas des conductivités du milieu. En revanche, si $V$ est solution de l'équation de Fredholm avec des valeurs de conductivité $\sigma_k$, $k \in [1..nbsurf]$, alors $\frac{1}{\lambda}V$ est également solution de l'équation de Fredholm lorsque les conductivités des milieux ont pour valeurs $\lambda\sigma_k$, $k \in [1..nbsurf]$. Le potentiel électrique est donc inversement proportionnel à un facteur multiplicatif affectant toutes les conductivités.

<footnote n="0">⁰Cette annexe est référencée aux sections 1.2.2, page 16, et 4.1, page 94</footnote>
## 2 Effet d'un facteur multiplicatif affectant toutes les conductivités en EEG

Le champ magnétique est donné par la loi de Biot et Savard (1.2.1.6 page 14) :

$$ \vec{B}(\vec{r}) = \vec{B}_o(\vec{r}) - \frac{\mu_o}{4\pi} \int_G \sigma(\vec{r}') \vec{\nabla} V \wedge \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} dv' $$

où

$$ \vec{B}_o(\vec{r}) = \frac{\mu_o}{4\pi} \int_G \vec{J}_p(\vec{r}') \wedge \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} dv' $$

et V est solution de l'équation de Fredholm.

$\vec{B}_o$ ne dépend pas des conductivités du milieu. Si V est solution de l'équation de Fredholm avec des valeurs de conductivité $\sigma_k$, $k \in [1..nbsurf]$, et $\vec{B}$ le champ magnétique donné par la loi de Biot et Savard pour cette valeur de V, alors lorsque les conductivités des milieux sont $\lambda\sigma_k$, $k \in [1..nbsurf]$, la solution de l'équation de Fredholm est $\frac{1}{\lambda}V$ et le champ magnétique donné par la loi de Biot et Savard est à nouveau $\vec{B}$. Le champ magnétique n'est pas affecté par un facteur multiplicatif de toutes les conductivités.
# Annexe E

## Forme matricielle de l'approche problème isolé

### 1 Introduction

L'approche problème isolé consiste à résoudre dans un premier temps un problème simplifié où le modèle est constitué de $p<n$ couches, puis d'utiliser la solution ainsi trouvée pour modifier le second membre du problème initial à n couches, afin d'obtenir une meilleure précision numérique. Sous la forme matricielle, l'équation du problème à n couches s'écrit :

$$A.V = V_o$$

Pour chaque couche, on peut définir les blocs correspondants des matrices $A$ et $V_o$, on obtient alors le découpage suivant :

$$A = \begin{pmatrix} A_{11} & \dots & A_{1p} & \dots & A_{1n} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ A_{p1} & \dots & A_{pp} & \dots & A_{pn} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ A_{n1} & \dots & A_{np} & \dots & A_{nn} \end{pmatrix} \quad \text{et} \quad V_o = \begin{pmatrix} V_{o1} \\ \vdots \\ V_{op} \\ \vdots \\ V_{on} \end{pmatrix}$$

Pour le problème à p couches, les conductivités seront notées $\tilde{\sigma}$ et pour le problème à n couches, $\sigma$. On a alors :

$$\tilde{\sigma}_i^+ = \tilde{\sigma}_i^- = 0 \text{ pour } i = 1, \dots, p-1$$

$$\tilde{\sigma}_p^+ = 0, \tilde{\sigma}_p^- = \sigma_p^-, \text{ et}$$

$$\tilde{\sigma}_i^+ = \sigma_i^+ \text{ et } \tilde{\sigma}_i^- = \sigma_i^- \text{ pour } i = p+1, \dots, n.$$

<footnote n="0">⁰Cette annexe est référencée à la section 1.2.3.4, page 20</footnote>
## 2 Problème à p couches.

Le potentiel électrique est régi en tout point par l'équation :

$$ \frac{\tilde{\sigma}_i^+ \Omega_{i,\vec{r}}^+ + \tilde{\sigma}_i^- \Omega_{i,\vec{r}}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} V_1(\vec{r}) = \frac{4\pi}{2\pi(\sigma_i^+ + \sigma_i^-)} V_o(\vec{r}) - \sum_{k=p}^{n} \frac{\tilde{\sigma}_k^- - \tilde{\sigma}_k^+}{2\pi(\sigma_i^+ + \sigma_i^-)} \int_{S_k \setminus \exists \vec{r}} V_1(\vec{r}') \vec{\nabla} \frac{1}{R} d\vec{S}'_k $$

qui se met sous la forme matricielle :

$$ MU_1 \cdot V_1 = V_{o1} $$

où :

$$ MU_1 = \begin{pmatrix} \frac{\sigma_p^+ + \sigma_p^-}{\sigma_p^- - \sigma_p^+} A_{pp} & A_{p,p+1} & \cdots & A_{pn} \\ A_{p+1,p} & A_{p+1,p+1} & & \vdots \\ \vdots & & \ddots & \vdots \\ A_{np} & \cdots & \cdots & A_{nn} \end{pmatrix} + D \text{ et } V_{o1} = \begin{pmatrix} \frac{\sigma_p^- + \sigma_p^+}{\sigma_p^-} V_{o_p} \\ V_{o_{p+1}} \\ \vdots \\ V_{o_n} \end{pmatrix} $$

D est une matrice diagonale donnée par $\sum_{j=p}^{n} (MU_1)_{ij} = 0, i = p \dots n$

## 3 Passage au problème à n couches

Pour le problème à n couches on pose :

$$ V(\vec{r}) = V_1(\vec{r}) + V_2(\vec{r}) $$

Où $V_1$ nous est donné par le problème isolé, et l'on cherchera donc à calculer $V_2$. Écrivons les équations régissant $V$ en remplaçant le terme $V_o$ par son expression dans le problème isolé.

$$ \frac{\sigma_i^+ \Omega_{i,\vec{r}}^+ + \sigma_i^- \Omega_{i,\vec{r}}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} V(\vec{r}) = \frac{4\pi}{2\pi(\sigma_i^+ + \sigma_i^-)} V_o(\vec{r}) - \sum_{k=1}^{n} \frac{\sigma_k^- - \sigma_k^+}{2\pi(\sigma_i^+ + \sigma_i^-)} \int_{S_k \setminus \exists \vec{r}} V(\vec{r}') \vec{\nabla} \frac{1}{R} d\vec{S}'_k $$

avec :

$$ \frac{4\pi}{2\pi(\sigma_i^+ + \sigma_i^-)} V_o(\vec{r}) = \frac{\tilde{\sigma}_i^+ \Omega_{i,\vec{r}}^+ + \tilde{\sigma}_i^- \Omega_{i,\vec{r}}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} V_1(\vec{r}) + \sum_{k=p}^{n} \frac{\tilde{\sigma}_k^- - \tilde{\sigma}_k^+}{2\pi(\sigma_i^+ + \sigma_i^-)} \int_{S_k \setminus \exists \vec{r}} V_1(\vec{r}') \vec{\nabla} \frac{1}{R} d\vec{S}'_k $$

d'où :

$$ \frac{\sigma_i^+ \Omega_{i,\vec{r}}^+ + \sigma_i^- \Omega_{i,\vec{r}}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} V_2(\vec{r}) = \left[ \begin{aligned} & \frac{\tilde{\sigma}_i^+ \Omega_{i,\vec{r}}^+ + \tilde{\sigma}_i^- \Omega_{i,\vec{r}}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} V_1(\vec{r}) - \frac{\sigma_i^+ \Omega_{i,\vec{r}}^+ + \sigma_i^- \Omega_{i,\vec{r}}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} V_1(\vec{r}) + \\ & \sum_{k=p}^{n} \frac{\tilde{\sigma}_k^- - \tilde{\sigma}_k^+}{2\pi(\sigma_i^+ + \sigma_i^-)} \int_{S_k \setminus \exists \vec{r}} V_1(\vec{r}') \vec{\nabla} \frac{1}{R} d\vec{S}'_k - \sum_{k=1}^{n} \frac{\sigma_k^- - \sigma_k^+}{2\pi(\sigma_i^+ + \sigma_i^-)} \int_{S_k \setminus \exists \vec{r}} V_1(\vec{r}') \vec{\nabla} \frac{1}{R} d\vec{S}'_k \\ & - \sum_{k=1}^{n} \frac{\sigma_k^- - \sigma_k^+}{2\pi(\sigma_i^+ + \sigma_i^-)} \int_{S_k \setminus \exists \vec{r}} V_2(\vec{r}') \vec{\nabla} \frac{1}{R} d\vec{S}'_k \end{aligned} \right] $$
On appellera par la suite les termes entre crochets $V_{o2}(\vec{r})$. On a alors :

$$V_{o2}(\vec{r}) = -\frac{\delta_{i,p} \sigma_p^+ \Omega_{p,\vec{r}}^+}{2\pi(\sigma_i^+ + \sigma_i^-)} V_1(\vec{r}) + \frac{\sigma_p^+}{2\pi(\sigma_i^+ + \sigma_i^-)} \int_{S_p \setminus \exists \vec{r}} V_1(\vec{r}') \vec{\nabla} \frac{1}{R} \cdot d\vec{S}_p'$$

De plus, si on applique A à :

$$\begin{cases} \tilde{V}_1(\vec{r}) = V_1(\vec{r}) & \text{si } \vec{r} \in S_p \\ \tilde{V}_1(\vec{r}) = 0 & \text{sinon} \end{cases}$$

on obtient :

$$A \tilde{V}_1(\vec{r}) = \frac{\sigma_p^- - \sigma_p^+}{2\pi(\sigma_i^+ + \sigma_i^-)} \int_{S_p \setminus \exists \vec{r}} \tilde{V}_1(\vec{r}') \vec{\nabla} \frac{1}{R} \cdot d\vec{S}_p' + \delta_{i,p} \frac{\sigma_p^+ \Omega_{p,\vec{r}}^+ + \sigma_p^- \Omega_{p,\vec{r}}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} \tilde{V}_1(\vec{r})$$

On peut alors écrire :

$$V_{o2}(\vec{r}) = \frac{\sigma_p^+}{\sigma_p^- - \sigma_p^+} A \tilde{V}_1(\vec{r}) - \delta_{i,p} \frac{\sigma_p^+}{\sigma_p^- - \sigma_p^+} \frac{\sigma_p^+ \Omega_{p,\vec{r}}^+ + \sigma_p^- \Omega_{p,\vec{r}}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} \tilde{V}_1(\vec{r}) - \frac{\delta_{i,p} \sigma_p^+ \Omega_{p,\vec{r}}^+}{2\pi(\sigma_i^+ + \sigma_i^-)} V_1(\vec{r})$$

$$V_{o2}(\vec{r}) = \frac{\sigma_p^+}{\sigma_p^- - \sigma_p^+} A \tilde{V}_1(\vec{r}) - \delta_{i,p} \frac{\sigma_p^+}{\sigma_p^- - \sigma_p^+} \left( \frac{\sigma_p^+ \Omega_{p,\vec{r}}^+ + \sigma_p^- \Omega_{p,\vec{r}}^- + (\sigma_p - \sigma_p^+) \Omega_{p,\vec{r}}^+}{2\pi(\sigma_i^+ + \sigma_i^-)} \right) V_1(\vec{r})$$

$$V_{o2}(\vec{r}) = \frac{\sigma_p^+}{\sigma_p^- - \sigma_p^+} A \tilde{V}_1(\vec{r}) - \delta_{i,p} \frac{\sigma_p^+}{\sigma_p^- - \sigma_p^+} \left( \frac{\sigma_p^- (\Omega_{p,\vec{r}}^- + \Omega_{p^+,\vec{r}}^-) + \sigma_p^+ (\Omega_{p,\vec{r}}^+ - \Omega_{p^+,\vec{r}}^+)}{2\pi(\sigma_i^+ + \sigma_i^-)} \right) V_1(\vec{r})$$

$$V_{o2}(\vec{r}) = \frac{\sigma_p^+}{\sigma_p^- - \sigma_p^+} A \tilde{V}_1(\vec{r}) - \frac{2\sigma_p^+ \sigma_p^-}{(\sigma_p^- - \sigma_p^+)(\sigma_p^+ + \sigma_p^-)} \tilde{V}_1(\vec{r})$$

$V_2$ est alors solution de l'équation matricielle :

$$AV_2 = V_{o2}$$

avec

$$V_{o2} = \frac{\sigma_p^+}{\sigma_p^- - \sigma_p^+} \begin{pmatrix} A_{11} & \dots & A_{1p} & \dots & A_{1n} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ A_{p1} & \dots & A_{pp} & \dots & A_{pn} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ A_{n1} & \dots & A_{np} & \dots & A_{nn} \end{pmatrix} \begin{pmatrix} 0 \\ 0 \\ V_{1p} \\ 0 \\ 0 \end{pmatrix} - \frac{2\sigma_p^+ \sigma_p^-}{\sigma_p^- - \sigma_p^+} \begin{pmatrix} 0 \\ 0 \\ V_{1p} \\ 0 \\ 0 \end{pmatrix}$$
<figure><img src="image_79.png" /></figure>
# Annexe F

## Calcul des matrices dans l'approche LINT
BEM avec des triangles courbes

### 1 Approche LINT sur un triangle plat

La matrice s'écrit :

$$A_{m,n} = \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} \delta_{m,n} + \frac{1}{2\pi} \frac{\sigma_j^- - \sigma_j^+}{\sigma_i^+ + \sigma_i^-} \sum_{\vec{r}_n \in \Delta_{kln}} \int_{\Delta_{kln} - \exists \vec{r}_m} \phi_n(\vec{r}) \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} d\vec{S}$$

et l'intégrale est calculée analytiquement avec les formules de De Munck qui sont données en annexe C.

### 2 Approche LINT sur un triangle courbe

Le triangle courbe K est défini comme suit : soit $\hat{K}$ le triangle plat de référence qui a pour sommets (0,0), (0,1) et (1,0). Il existe une transformation $F_{\hat{K},K}$ de l'espace affine, polynômiale de degré 2, telle que $F_{\hat{K},K}(\hat{K}) = K$.

<figure><img src="image_80.png" /></figure>

<footnote n="0">⁰Cette annexe est référencée à la section 2.2, page 54</footnote>
On peut alors écrire l'équation paramétrique du triangle dans $\mathbb{R}^3$ sous la forme :

$$F_{\hat{K},K}(u,v) = f_1(u,v)x_1 + f_2(u,v)x_2 + f_3(u,v)x_3 + f_{12}(u,v)x_{12} + f_{13}(u,v)x_{13} + f_{23}(u,v)x_{23}$$

avec :

$$
\left\{
\begin{array}{ll}
0 \le u + v \le 1, \\
u \in [0 \ 1], \\
v \in [0 \ 1], \\
K \quad & \text{a pour sommets } x_1, x_2, x_3, \\
K \quad & \text{a pour milieux de ses côtés } x_{12}, x_{13}, x_{23}, \\
f_1(u, v) = (1 - u - v)(1 - 2u - 2v), \\
f_2(u, v) = u(2u - 1), \\
f_3(u, v) = v(2v - 1), \\
f_{12}(u, v) = 4u(1 - u - v), \\
f_{13}(u, v) = 4v(1 - u - v), \\
f_{23}(u, v) = 4uv.
\end{array}
\right.
$$

Les $f$ sont les fonctions de base P2 des éléments finis. Le calcul de l'intégrale sur $K$ est alors modifié, car pour une surface courbe les formules de De Munck ne s'appliquent pas. On propose alors de subdiviser le triangle $K$ en petits triangles plats sur lesquels les formules de De Munck s'appliquent.

<figure><img src="image_81.png" /></figure>

Le calcul de la matrice s'écrit :

$$A_{m,n} = \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} \delta_{m,n} + \frac{1}{2\pi} \frac{\sigma_j^- - \sigma_j^+}{\sigma_i^+ + \sigma_i^-} \sum_{\vec{r}_n \in \Delta_{kln}} \int_{\Delta_{kln} - \exists \vec{r}_m} \tilde{\phi}_n(\vec{r}) \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} \cdot d\vec{S}$$

Où

$$\tilde{\phi}_n(\vec{r}) = \phi_3(F_{\hat{K},\Delta_{kln}}^{-1}(\vec{r})), \forall \vec{r} \in \Delta_{kln}$$

Découpons, comme le propose De Munck, le triangle courbe $\Delta_{kln}$ en petits triangles plats. On a alors :

$$\int_{\Delta_{kln}} \tilde{\phi}_n(\vec{r}) \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} \cdot d\vec{S}' \approx \sum_{\Delta_{pqr} \subset \Delta_{kln}} \int_{\Delta_{pqr}} \tilde{\phi}_n(\vec{r}) \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} \cdot d\vec{S}'$$
avec

$$ \tilde{\phi}_n(\vec{r}) = \tilde{\phi}_n(\vec{r}_p)\phi_p(\vec{r}) + \tilde{\phi}_n(\vec{r}_q)\phi_q(\vec{r}) + \tilde{\phi}_n(\vec{r}_r)\phi_r(\vec{r}), \quad \forall \vec{r} \in \Delta_{pqr} $$

Ainsi, en stockant les $\tilde{\phi}_n(\vec{r}_i)$ lors de la construction du maillage, on peut calculer la matrice de la méthode des intégrales de surface avec l'équation suivante :

$$ A_{m,n} = \frac{\sigma_i^+ \Omega_{i,r_m}^+ + \sigma_i^- \Omega_{i,r_m}^-}{2\pi(\sigma_i^+ + \sigma_i^-)} \delta_{m,n} + \frac{1}{2\pi} \frac{\sigma_j^- - \sigma_j^+}{\sigma_i^+ + \sigma_i^-} \sum_{\vec{r}_n \in \Delta_{kln}} \sum_{\Delta_{pqr} \subset \Delta_{kln}} \sum_{i=p,q,r} \tilde{\phi}_n(\vec{r}_i) \int_{\Delta_{pqr}} \phi_i(\vec{r}) \vec{\nabla} \frac{1}{|\vec{r}_m - \vec{r}|} d\vec{S} $$

Cette fois, les intégrales peuvent bien être calculées avec les formules de De Munck.
<figure><img src="image_82.png" /></figure>
# Annexe G

## Article paru dans : Clinical

## Neurophysiology 110 (1999) 2176-2188

<footnote n="0">⁰Cette annexe est référencée à la section 3.1.2, page 62</footnote>
# An evaluation of dipole reconstruction accuracy with spherical and realistic head models in MEG.

A. Crouzeix<fnref n="1" />, B. Yvert, O. Bertrand, J. Pernier

Mental Processes and Brain Activation, INSERM U280, Lyon, France

## Keywords

Magnetoencephalography, Realistic head model, Spherical head model, Dipole localization, Boundary Element Method, inverse problem.

## Abstract

MEG forward problem has been solved for about 2000 dipoles placed on the brain surface using a very fine three-layer realistic model of the head and the Boundary Element Method (BEM). For each dipole, spherical models, one-layer realistic BEM models and coarser three-layer realistic BEM models were used to reconstruct the dipole. It was found that the localization bias induced by using a spherical model of the head increased from 2.5 mm in the upper part of the head to 12 mm in the lower part, on average. It was also found that, for the same computing time, a three-layer model of the head gave on average 2 mm better localization errors than a one-layer model of the head. Orientation errors of less than 20° could only be retrieved with a three-layer realistic model. Localization and orientation errors highly depended on the dipole position in the brain.

## 1 Introduction

The spherical head model is widely used for estimating sources of the electrical brain activity from magnetoencephalographic signals (MEG). The general belief that spherical model is more accurate for MEG than for electroencephalographic signals (EEG) stems from the fact that MEG does not depend on conductivity ratios between layers in a spherical geometry. The bias induced on dipole reconstruction by the spherical approximation has been systematically tested for magnetic data by Fuchs et al. [1998], but only using a 1-layer realistic model of the head.

More sophisticated spherical models consisting of eccentric spheres have been proposed in MEG by [Meijs and Peters, 1987a], and further studied by [Cuffin, 1991] in terms of varying

<footnote n="1">¹Corresponding author. Anne Crouzeix, INSERM - Unité 280, 151, Cours Albert Thomas, F-69424 Lyon cedex 03, France. Tel +33-(0)4 72 68 19 01; fax +33-(0)4 72 68 19 02. E-mail: crouzeix@lyon151.inserm.fr</footnote>
skull and scalp thicknesses. The difference in localization accuracy between the concentric sphere model and these models was found to be negligible (0.3 mm error). However, only one tangential dipole placed on an axis of symmetry of the model geometry was considered.

More realistic head models have been introduced by Hämäläinen [1987]; Meijs et al. [1987b]. Magnetic resonance images (MRI) segmentation is used to obtain meshes for real-shaped skull, scalp and brain, and the Boundary Element Method (BEM) can then be applied to solve the forward problem. Since magnetic field is less sensitive to conductivity distribution, it is often considered that a model with only the brain layer is sufficient. In Zanow [1997], the localization error obtained on data generated by a one-layer realistic model with a spherical inverse model was ranging from 1 mm for superficial dipoles near the vertex to 25 mm for deeper frontal dipoles. In Tomita et al. [1996], a similar simulation gave localization errors ranging from 0.67 mm for a shallow occipital dipole to 21 mm for the deepest frontal dipole. The orientation error for almost tangential dipoles was ranging from 1° to 63°. However, 12 of the 15 dipoles had orientation errors less than 10°. In Fuchs et al. [1998], a systematic simulation study with about 3600 dipoles and a 1-layer realistic model for generating the data showed that the localization error using a spherical model ranged from 0.4 mm to 45.1 mm for tangential dipoles in a 122 channels helmet system. Smaller localization errors where obtained for tangential dipoles located near the pick-up coils.

Beside those simulation studies, the spherical model accuracy was also tested experimentally in MEG, using realistic phantoms of the head, or implanted dipoles in epileptic patients undergoing presurgical intracerebral recordings.

In Barth et al. [1986], a human cadaver head in which the brain had been replaced with salt jelly was used. Five dipoles and 81 sequentially recorded positions for the single MEG sensor were used. Localization errors obtained with a spherical model ranged from 1.3 mm for shallow sources to 9.5 mm for the deepest sources. In Menninghaus et al. [1994b], a nonconducting plastic skull filled with a saline solution was used, corresponding to a one-layer physical model of the head. Localization errors ranging from 3 mm to 9 mm were found for 4 dipoles with the 37 sensor configuration of the BTI system (MAGNES(R)). In Leahy et al. [1998], a human skull was impregnated and filled with saline gelatin, and coated with latex to obtain a three-layer physical model of the head. Magnetic fields successively generated by 32 dipoles (left central sulcus, left calcarine fissure, and frontal area) were recorded with a Neuromag-122 planar gradiometer. Using a spherical model, and the R-MUSIC localization method, they found errors ranging from 1 mm to 17 mm. The BEM, based on realistically deformed spherical meshes, lead to a 0.5-mm improvement of the localization error only.

In Cohen et al. [1990]; Cohen and Cuffin [1991]; Balish et al. [1991] MEG localization errors using implanted dipoles were found to be ranging from 3 mm to 20 mm, with a mean
varying from 8 to 10 mm, strongly depending on the dipole depth, and provided that the signal-to-noise ratio (SNR = $\sqrt{\frac{\langle B(t)^2 \rangle_{\text{dipole on}}}{\langle B(t)^2 \rangle_{\text{dipole off}}}}$) was greater than 10. However, SNR was acceptable for only a few dipoles, possibly because most dipoles had a radial orientation with respect to the inner skull surface. InRose et al. [1991] 3 implanted dipoles in the infero-temporal lobe were localized within 1-4 mm.

The localization errors found in these simulation and experimental studies using the spherical head model for dipole reconstruction are summarized in Table G.1. However, most of these MEG studies considered very few dipole locations.

Yvert et al. [1997] reported a systematic evaluation of the dipole localization accuracy obtained from electrical signals with a spherical model. The work presented here aims at extending this study to the magnetic case.

For a given 143 MEG sensor configuration, magnetic fields were computed with a very fine realistic 3-layer BEM model for 2003 dipoles distributed on the cortical surface. Then, the localization bias due to the use of a spherical head model, a one-layer realistic BEM model, or a coarse three-layer realistic BEM model were studied. Dipole orientation errors were also systematically evaluated in the case of realistic geometries. This aspect was less frequently reported in the literature [Tomita et al., 1996; Leahy et al., 1998].

This MEG simulation study focused on the evaluation of the localization bias due to different level of geometrical approximations of the fine reference model. The influence of other sources of errors such as background recording noise or coil location uncertainties are beyond the scope of this paper.
<table><caption>TAB. G.1: Localization errors found in the literature using the concentric sphere model for inverse computations and simulated or real magnetic data. SNR : signal to noise ratio.</caption><thead><tr><th>Authors</th><th>MEG data</th><th>Number of dipoles</th><th>Dipole depth (mm)</th><th>Dipole orientation</th><th>Localization error (mm)</th></tr></thead><tbody><tr><td>Cuffin [1991]</td><td>simulated data, ecc. sph. model</td><td>1</td><td>31</td><td>Tangential</td><td>0.3</td></tr><tr><td>Zanow [1997]</td><td>simulated data, one-layer realistic model</td><td>30</td><td>7 - 53</td><td>Tangential</td><td>1 - 25</td></tr><tr><td>Tomita et al. [1996]</td><td>simulated data, one-layer realistic model</td><td>15</td><td>≈ 15-55</td><td>Tangential</td><td>0.7 - 21.4</td></tr><tr><td>Fuchs et al. [1998]</td><td>simulated data, one-layer realistic model</td><td>≈ 3600</td><td>≈ 0-80</td><td>1/3 Tangential, 1/3 intermediate, 1/3 radial</td><td>0.4-45.1, 0.2-70.9, 13.1-87.9</td></tr><tr><td>Barth et al. [1986]</td><td>Human cadaver</td><td>5</td><td>≈ 20 - 75</td><td>Variable</td><td>1.3 - 9.5</td></tr><tr><td>Menninghaus et al. [1994b]</td><td>Phantom, one-layer</td><td>4</td><td>10 - 30</td><td>Tangential</td><td>3 - 9</td></tr><tr><td>Leahy et al. [1998]</td><td>Phantom, 3-layer</td><td>32</td><td>10 - 41</td><td>Variable</td><td>1 - 17 (SNR &lt; 4)</td></tr><tr><td>Balish et al. [1991]</td><td>Implanted dipoles</td><td>19</td><td>20 - 63</td><td>Tangential</td><td>3 - 40 3-20 for SNR &gt; 10</td></tr><tr><td>Cohen et al. [1990]</td><td>Implanted dipoles</td><td>12</td><td>≈ 30 - 45</td><td>10 radial 2 tangential</td><td>4 - 40 4-13 for SNR &gt; 1</td></tr><tr><td>Rose et al. [1991]</td><td>Implanted dipoles</td><td>3</td><td>≈ 15 - 45</td><td>Tangential</td><td>1 - 4</td></tr></tbody></table>
# 2 Methods

## 2.1 Mesh generation

Meshes for the BEM were derived from a stack of MRI slices following a three-step method : First, from each MRI slice, scalp, skull and brain contours were extracted. Second, contours were meshed together following the method described in [Yvert et al., 1995] leading to very refined meshes (about 20000 nodes). Third, mesh density was reduced according to the simplification technique proposed by Algorri and Schmitt [1996] in order to obtain a tolerable number of triangles for practical use (Table G.2).

## 2.2 Forward model

The simulated MEG data were computed using the BEM with linear interpolation of the potential on the surface elements [De Munck, 1992; Ferguson et al., 1994] and the isolated problem approach [Meijs et al., 1987b; Hämäläinen, 1987]. The model used for this computation (called REF model in the following) had about 4500 nodes and a high triangle density (mean : 3 triangles/cm²)( Table G.2) and is shown in Fig. G.1. The normalized conductivity values were 1,1/80,1 for the scalp, skull and brain compartments respectively, as suggested in Rush and Driscoll [1968] and Homma et al. [1995].

<table><caption>TAB. G.2: Characteristics of the models used for the forward and inverse problem.</caption><thead><tr><th>Model name</th><th>Model type</th><th>Number of layers</th><th>Number of nodes<br>per layer</th><th>Triangle density<br>(tri./cm²)<br>scalp-skull-brain</th></tr></thead><tbody><tr><td colspan="5">Forward problem</td></tr><tr><td>REF</td><td>Realistic</td><td>3</td><td>1500</td><td>2.1 - 3 - 4</td></tr><tr><td colspan="5">Inverse problem</td></tr><tr><td>GLSPH</td><td>Global<br>best-fitting<br>sphere</td><td>-</td><td>-</td><td>-</td></tr><tr><td>LOCSPH</td><td>Local best-fitting<br>sphere</td><td>-</td><td>-</td><td>-</td></tr><tr><td>1L900</td><td>Realistic</td><td>1</td><td>900</td><td>2.6</td></tr><tr><td>3L300</td><td>Realistic</td><td>3</td><td>300</td><td>0.5 - 0.7 - 1.2</td></tr><tr><td>3L900</td><td>Realistic</td><td>3</td><td>900</td><td>1.9 - 1.8 - 2.6</td></tr></tbody></table>
<figure><img src="image_83.png" /><figcaption>FIG. G.1: Simulation procedure. Models used for the forward and inverse procedure are presented with superimposed pickup coils positions. REF is a three-layer realistic model with a triangle density of 3 tri/cm², and is used for creating the simulated MEG data. GLSPH and LOCSPH are spherical models used for inverse computations. 1L900 is a one-layer realistic model with 2.6tri/cm², 3L300 and 3L900 are three-layer models with 0.7tri/cm² and 1.8tri/cm² respectively.</figcaption></figure>
A linear interpolation method being used for the BEM, a key feature of the meshes considered in this study is the number of nodes which determines directly the size of the matrix equations to be solved.

## 2.3 MEG sensor configuration

MEG data were computed on a whole-head sensor configuration from the CTF company (143 first order gradiometers, with a baseline coil 5 cm away from the pickup coils). This arrangement was realistically positioned around the REF model (Fig. G.1). The distance between scalp and sensors was as much constant as possible (minimum distance : 19 mm).

## 2.4 Choice of the dipole sites

The 2003 dipole sites were identical as those used in Yvert et al. [1997], enabling further comparison between EEG and MEG results. They were randomly chosen from the intersection of a cortical surface obtained by MRI segmentation and 10 horizontal planes separated by 10 mm (Fig. G.2). The distance between two neighboring dipoles within a plane was about 5 mm. Each dipole orientation was normal to the cortical surface at the dipole location. The dipole depth under the cortical surface was ranging from 0.9 to 46.9 mm (mean 16.4 mm ± 9.8 mm). As shown in Fig. G.2, deeper dipoles lead to magnetic fields having a lower root mean square value (RMS).

To characterize the dipole orientation, a "tangentiality" angle (between 0° and 90°) was defined as the angle between the dipole moment and its radial component in the global sphere model (GLSPH, described in the next section). In the presentation of the results, dipoles having a tangentiality angle between 0 and 5° will be referred to as radial dipoles, whereas those with a tangentiality angle between 85 and 90° will be referred to as tangential dipoles. Of the 2003 dipoles, 63% have a "tangentiality" angle between 60° and 90°. Figure G.2 gives an overview of the tangentiality over all 2003 dipoles.
<figure><img src="image_84.png" /></figure>

RMS

<figure><img src="image_85.png" /></figure>

Tangentiality

<figure><img src="image_86.png" /><figcaption>FIG. G.2: Top : A 3D view of the segmented cortical surface has been superimposed on a mid-sagital MRI slice. The projection of the dipoles positions on the segmented cortical surface are represented by red dots. Dipoles are placed on 10 horizontal slices spaced by 10 mm. Middle : Root Mean Square values (RMS) of the magnetic field created by each dipole. The dipole population is divided into four equally numbered parts depending on their associated RMS values. For each dipole position, a color dot is superimposed on MRI slices according to the RMS value for that dipole. Red : 0 fT<RMS<0.5 fT ; yellow : 0.5 fT<RMS<0.9 fT ; green : 0.9 fT<RMS<1.3 fT ; blue :1.3 fT<RMS<14 fT. Bottom : Tangentiality values. For each dipole position, a color dot is superimposed on MRI slices, the color being chosen according to the tangentiality of that dipole : red : 0-30 degrees (almost radial) ; yellow : 30-60 degrees ; green : 60-90 degrees (almost tangential).</figcaption></figure>
## 2.5 Inverse models

### 2.5.1 The spherical models

Four spherical models were considered for solving the inverse problem :

1. A global sphere model (GLSPH, Fig. G.1) which best fitted the surface of the head-shape covered by all coils.

2. A brain-fitted global sphere model (BRAINGLSPH), similar to GLSPH but fitting the brain surface instead of the head shape.

3. A local sphere model (LOCSPH, Fig. G.1) which best fitted a surface of the head-shape near the coils where the signal had its minima and maxima. This surface was determined as follows :

(a) The distance *d* between the pickup coils recording the minimum and maximum signals was computed.

(b) Then, a subset of coils whose distance from the middle of those two coils was less than $2 * d$ was selected. If more than half of the total number of coils were selected, only the global sphere was considered, and the dipole was discarded from the comparison between spherical inverse models.

(c) The local sphere best fitted the surface of the head-shape whose distance from the middle of those two coils was less than $2 * d$.

4. A brain-fitted local sphere model (BRAINLOCSPH), similar to LOCSPH but fitting the brain surface instead of the head shape.

This procedure for selecting a subset of coils and for defining a local sphere has been chosen because it could usually be applied in practice, when a head-shape discretization is available : it does not require any a priori assumption on the possible location of the sources. Spheres fitting the head-shape are usually chosen only when no MRI is available. In the other case, brain fitted spheres are usually preferred.

Apart from those four spherical models, the sensitivity of spherical models to sphere center position was tested using six sphere models shifted from BRAINLOCSPH by 1 cm in the six axial directions.

### 2.5.2 The one-layer realistic model

From a topographical analysis of forward problems solved with one-layer and three-layer models for a few source configurations, Hämäläinen and Sarvas [1989] have suggested that since the skull has a very high resistivity compared to scalp and brain tissues, i.e. major currents are flowing in the brain, only one layer is needed to model the head. Here, we
extended this study by systematically testing this assumption for about 2000 dipoles. We used a one-layer model (Fig. G.1, Table G.2.) having about 900 nodes and a density of 2.5 triangles/cm².

### 2.5.3 The three-layer realistic model

We also considered two realistic three-layer models : one coarse (3L300), and one fine (3L900) having respectively 300 and 900 nodes per layer (see Fig. G.1, Table G.2). The most time-consuming steps in the inverse computations being matrix inversion and matrix product, 3L300 and 1L900 had similar computational time.

## 2.6 The inverse procedure

The inverse problem was solved iteratively by varying the dipole position and orientation parameters p so as to minimize the following criterion :

$$ \epsilon = \sum_{i=1}^{N} (B_{REF}^{i} - B_{INVmodel}^{i}(\mathbf{p}))^{2} $$

where N is the number of channels, and $B_{REF}^{i}$ is the magnetic field computed with the REF model at the $i^{th}$ channel, and $B_{INVmodel}^{i}(\mathbf{p})$ is the magnetic field produced at the $i^{th}$ channel by a dipole with parameters p for a given inverse model.

A mixed linear/non-linear iterative algorithm (Marquardt) was used for the minimization [Scherg, 1990; Verkindt et al., 1995], which required initial conditions for the dipole parameters. In this study, the exact position of the original dipole was taken as the initial conditions. For 186 shallow dipoles, the 3L300 mesh was so rough that these dipoles fell outside the meshed brain surface, and were thus discarded from the statistical analysis. Beside this, for 50 dipoles, the relative residual error ($\epsilon$ divided by the mean square value of the signal) at the end of the inverse procedure was over 30%, indicating that no reasonable solution was found. These cases, which mostly happened with spherical models, were also discarded.

## 2.7 Error criteria

The localization error $D$ is defined as the distance in millimeters between the reconstructed dipole and the original dipole, the orientation error $\theta$ as the angle between the reconstructed dipole moment and the original dipole moment, and the moment amplitude error $\Delta M$ as the difference between the reconstructed dipole moment amplitude and the original dipole moment amplitude. $\Delta M$ is expressed in percentage of the original dipole
moment amplitude. Negative and positive values indicate underestimation and overestimation, respectively of the dipole strength.
<figure><img src="image_87.png" /><figcaption>FIG. G.3: Mean localization and orientation errors for spherical models with respect to slice number (A) and tangentiality angle (B). Error bars represent the standard error of the mean. Only the 824 dipoles for which the local sphere is different from the global sphere and leading to a residual fitting error less than 30% are plotted.</figcaption></figure>
<figure><img src="image_88.png" /><figcaption>FIG. G.4: Mean localization (A), orientation (B), and moment amplitude (C) errors with respect to the tangentiality angle. Error bars represent the standard error of the mean. Only the 1768 dipoles inside all meshed brain layers and leading to a residual fitting error less than 30% are considered.</figcaption></figure>
<figure><img src="image_89.png" /><figcaption>FIG. G.5: Localization error for BRAINLOCSPH and realistic models. (Note that 1L900 and 3L300 have the same computing cost). For each dipole position, a color dot is superimposed on MRI slices, the color being chosen according to the localization error obtained for that dipole : red : 0-3 mm ; yellow : 3-6 mm ; green : 6-12 mm ; blue : >12 mm.</figcaption></figure>
# 3 Results

## 3.1 Bias due to the spherical approximation

For 824 dipoles (depth usually inferior to 2 cm), the local sphere (LOCSPH) and the global sphere (GLSPH) were different. For all these dipoles the brain fitted local sphere (BRAINLOCSPH) and the brain fitted global sphere (BRAINGLSPH) were also different. The distance between the centers of the scalp-fitting spheres ranged from 2.5 mm to 58 mm (mean 21 mm ± 11.5 mm), and for brain fitting spheres, from 2.7 to 39 mm (mean 14.7 mm ± 6.9 mm). For the comparison of GLSPH, BRAINGLSPH, LOCSPH and BRAINLOCSPH, only these 824 dipoles were considered.

### 3.1.1 Localization errors

We found that for those 824 dipoles, GLSPH, LOCSPH, BRAINGLSPH and BRAINLOCSPH lead to a mean localization error of 3.5 mm, 2.4 mm, 3.1 mm and 2.5 mm, respectively, in the upper part of the head (slices 1-5 in Fig. G.3 A and Fig. G.5). The 1.1 mm difference between GLSPH and LOCSPH was significant (paired t test: *p* < 0.0001), and so was the 0.6 mm difference between BRAINGLSPH and BRAINLOCSPH. The 0.1 mm difference between LOCSPH and BRAINLOCSPH was not significant (*p* = 0.29).

In the lower part of the head (slices 6-10), GLSPH, LOCSPH, BRAINGLSPH, and BRAINLOCSPH lead to increased mean localization errors of 5.3 mm, 8.7 mm, 4.6 mm and 4 mm, respectively. This vertical slice effect was significant (*p* < 0.0001). A more detailed analysis in this region showed that:

i. In the posterior part of the brain, all spherical models except BRAINLOCSPH gave a mean localization error between 3.7 mm and 3.9 mm: there was no statistical difference between these models (*p* > 0.6). BRAINLOCSPH gave a mean localization error of 2.9 mm which was statistically different from other models (*p* < 0.0001).

ii. In the anterior part of the brain, GLSPH, LOCSPH, BRAINGLSPH, and BRAINLOCSPH lead to a mean localization error of 6.5 mm, 13.4 mm, 5.2 mm and 5 mm respectively. All those differences were significant (*p* < 0.0001) except between BRAINLOCSPH and BRAINGLSPH (*p* = 0.62).

When varying the BRAINLOCSPH sphere center in the six axial directions, localization errors varied on average of 0.8 mm. No particular direction gave systematically better results.
<figure><img src="image_90.png" /><figcaption>FIG. G.6: Orientation error for spherical and realistic models. For each dipole position, a color dot is superimposed on MRI slices, the color being chosen according to the orientation error obtained for that dipole: a red dot indicates an error between 0 and 10 degrees, a yellow dot, between 10 and 20 degrees, a green dot, between 20 and 40 degrees, and a blue dot, over 40 degrees.</figcaption></figure>
<figure><img src="image_91.png" /><figcaption>FIG. G.7: Left column : Mean localization (A), orientation (B), and moment amplitude (C) errors with respect to slice number. Middle column : same for the 75% of the dipole population which has the greatest RMS. Right column : EEG results obtained by Yvert et al. (1997) are shown for comparison. In the latter study, 3L660 was a 3-layer realistic model with about 660 nodes per layer called IUMM. GLSPH fitted all electrode positions and was called Model 2. Error bars represent the standard error of the mean. Only the 1768 dipoles inside all meshed brain layers and leading to a residual error less than 30% are considered for MEG results.</figcaption></figure>
### 3.1.2 Orientation errors

Orientation errors were very similar for all spherical models (Fig. G.3 A, Fig. G.6), and no obvious slice effect could be seen. It should be noted that each slice included a combination of dipoles with various orientations. As expected, radial dipole leads to orientation errors of 90°, whereas this error fell to 10° for tangential dipoles (Fig. G.3 B).

## 3.2 Errors obtained with realistic models

All statistics included in this part were computed with the 1768 dipoles which fell inside the meshed brain layer for all meshes considered and gave a residual fitting error below 30%. Results were compared to those obtained with GLSPH and BRAINLOCSPH.

### 3.2.1 Localization errors

Localization errors obtained with the different realistic models (1L900, 3L300 and 3L900) are given in Fig. G.4 A, Fig. G.5, Fig. G.7 A

* For 1L900 (Fig. G.7 A), the mean localization errors increased from 2.4 mm in the upper part of the head (slices 1-5) to 4.9 mm in the lower part of the head (slices 6-10). They varied from 8.3 mm to 2.5 mm for radial and tangential dipoles respectively (Fig. G.4 A).

* For 3L300 (Fig. G.7 A) the mean localization errors increased from 1.3 mm in the upper and middle part of the head (slices 1-8), to 3.5 mm in the lower part of the head (slices 9-10). They varied between 4 mm and 0.97 mm for radial and tangential dipoles respectively (Fig. G.4 A).

* For 3L900 (Fig. G.7 A) the mean localization errors increased from 0.4 mm in the upper and central part of the head (slices 1-8), to 1.4 mm in the lower part of the head (slices 9-10). They varied from 2.6 mm to 0.4 mm for radial and tangential dipoles respectively (Fig. G.4 A).

All these differences between lower and upper parts of the head were statistically significant (t-test : *p* < 0.0001). Dipoles were localized better than 3 mm in 59%, 90%, and 98% of the cases with 1L900, 3L300, and 3L900, respectively (Fig. G.5).

### 3.2.2 Orientation errors

Orientation errors obtained with the different realistic models (1L900, 3L300 and 3L900) are given in Fig. G.4 B, Fig. G.6, Fig. G.7 B.
* For 1L900 the mean orientation errors remained between 13° and 31° when the dipole tangentiality varied.

* For 3L300 the mean orientation errors remained between 3° and 21° when the dipole tangentiality varied.

* For 3L900 the mean orientation errors remained between 2° and 8° when the dipole tangentiality varied.

Dipole orientations were retrieved with an accuracy better than 10° for 27%, 60%, and 87% of the dipoles with 1L900, 3L300, and 3L900, respectively (Fig. 6).

### 3.2.3 Moment amplitude error

Moment amplitude errors obtained with the different realistic models (1L900, 3L300 and 3L900) are given in Fig. G.4 C, Fig. G.7 C.

* For 1L900 the mean moment amplitude errors varied between -50% and 0% when the dipole tangentiality varied from 0 to 60 degrees, and between 0% and 24% when the dipole tangentiality varied from 70 to 90 degrees.

* For 3L300 the mean moment amplitude varied between -27% and 0% when the dipole tangentiality varied from 0 to 50 degrees, and between 0% and 6% when the dipole tangentiality varied from 60 to 90 degrees.

* For 3L900, the mean moment amplitude varied between -11% and 0% when the dipole tangentiality varied from 0 to 60 degrees, and between 0% and 2% when the dipole tangentiality varied from 70 to 90 degrees.

Dipole moment amplitudes were retrieved with an accuracy between -10% and 10% for 32%, 67%, and 87% of the dipoles with 1L900, 3L300, and 3L900, respectively.

## 4 Discussion

The goal of the present MEG simulation study was to estimate the localization bias induced by different approximations of a given realistic reference model (REF). In particular, this work aimed at determining the brain regions where a spherical model provides similar results to those obtained with realistic models, and at determining the influence of realistic model characteristics on source position, orientation, and strength reconstruction. It is very likely that implanted dipole studies (phantoms or epileptic patients) provide the best realistic data for evaluating source localization methods. However, for practical limitations, only a few tens of dipole positions are generally considered [Barth et al., 1986; Cohen et al., 1990; Balish et al., 1991; Rose et al., 1991; Menninghaus et al., 1994b; Leahy et al., 1998]. By contrast, in
simulation studies, like the present one, one has to choose a particular mathematical model of
the head to generate the data set. However, they offer a mean for systematically evaluating
localization errors all over the brain.

4.1 Localization bias due to spherical model

We have found that the localization bias induced by using a spherical model of the head ranged from about 2.5 mm in the uppermost part of the head up to about 8 mm in the temporal lobe. This "vertical position effect" has also been observed in EEG using simulated data [Yvert et al., 1997], which results are shown here for comparison in Fig. G.7. Similar important errors in the lower part of the head have been recently reported using spherical models for localization of equivalent sources underlying interictal spike discharges in epileptic patients [Merlet et al., 1998, Fig. 4]. This effect is much less pronounced with realistic models both in EEG and MEG. The spherical model induces a smaller localization bias in the most spherical-like regions of the head which are also well covered by whole head helmet systems. The range of localization errors is similar to what Zanow [1997] and Fuchs et al. [1998] found when using a 1-layer model of the head for the simulated forward computation, and agrees with the range of errors found with artificial dipoles in epileptic patients [Balish et al., 1991; Cohen et al., 1990; Rose et al., 1991] or phantoms [Barth et al., 1986; Menninghaus et al., 1994b]. At first sight, our results seem in contradiction with those of Leahy et al. [1998], who found that a very refined three-shell BEM model gave localization errors only 0.5 mm better than a local sphere, on average. However, their dipoles were distributed only along trajectories mimicking the central sulcus and upper occipital cortex, corresponding to brain areas where spherical and realistic models lead to similar errors in our simulation. Nevertheless, one of their dipoles was located in the lowest part of the head (dipole #12), for which the BEM model gave an error of about 4 mm, while the sphere model leads to about 17 mm error [Leahy et al., 1998, Fig. 11]. For the other dipoles, which were located upper in the head, realistic and spherical models lead to comparable errors. Hence, our simulations are consistent with their results. We also compared local and global sphere models. We found that for scalp-fitting spheres, the local sphere was slightly better (1 mm) than the global sphere in the upper part of the head. However, the situation was reversed in the lower part of the brain (14 mm for LOCSPH and 7 mm for GLSPH). By contrast, brain-fitted local spheres gave a significantly improved accuracy of 0.6 mm compared to brain-fitted global sphere. The local sphere has been systematically considered in the literature by authors who were using MEG devices with sensors partially covering the head-shape (BTI 37 channel system for instance). This approach was mainly used to explore the somatosensory and auditory cortices, in which areas we have found slightly smaller errors with the scalp-fitting local sphere. In more recent studies
[Mäkelä et al., 1994; Yvert et al., 1998] based on whole-head MEG systems, authors used a local sphere fitting a part of head shape covered by a subset of coils selected from a priori knowledge on the dipole location. Such choice might lead to important errors (over 8 mm) for sources in the infero-temporal or frontal regions, except if a brain-fitted local sphere is used. However, this sphere can only be chosen when MRIs are available, which is not often the case.

## 4.2 Influence of the number of layers with realistic models.

We have shown that the number of layers in BEM models may have an important influence on MEG single dipole reconstruction. We have found that the 3-shell model 3L300 leads to errors 2 mm smaller on average than with the 1-shell model 1L900 for the same computing time. Moreover, as discussed below, the 3-layer models allow a much better reconstruction of the dipole orientation and moment amplitude. This result may be mainly explained by a better account for secondary current sources with a 3-layer model. It should be recalled however that our data are not real but simulated with a 3-layer BEM model of the head. Hence, whether the very often proposed 1-layer realistic model is an oversimplification of the head geometry should be further investigated in more realistic situations such as implanted dipoles in epileptic patients.

## 4.3 Orientation and moment amplitude errors

We found that the dipole orientation and moment amplitude can be well estimated with realistic models of the head. It is commonly stated that MEG is not sensitive to the radial components of the dipoles. Our study showed that in realistic cases, the magnetic field distribution partially depends on this radial component, and that this information is sufficient (at least, in noiseless situations like here) to allow an estimation of the dipole orientation and moment amplitude. Fig. G.4 B shows that while the ability of the spherical model to reconstruct the dipole orientation and moment amplitude depends strongly on the dipole tangentiality (very poor for radial source), this is surprisingly not the case with realistic models. Indeed, for radial dipoles, the dipole orientation is reconstructed within less than 20° with both 3-layer realistic models, and 40° with the 1-layer realistic model. The moment amplitude is retrieved with less than 15 % error for radial dipoles and 3L900. For tangential dipoles the spherical models reconstruct the orientation of the dipoles within 10°. In this particular case it is even better than realistic models, most likely because the dipole is constrained to remain tangential, while the dipole can vary freely in orientation with realistic models. Fig. G.7 shows that 3-layer BEM models retrieve the dipole orientation within comparable errors between EEG and MEG. Although deep and radial dipoles are often
thought to be magnetically silent, in our data and others [Haueisen et al., 1995; Menninghaus and Lütkenhöner, 1994a], such sources still produce a RMS signal value being 20% of that produced by a tangential source. Hence, the ability of realistic model to retrieve radial and deep sources should be further investigated in noisy simulations, which was beyond the scope of this study. Nevertheless, we have here considered the error results obtained for only the 75% of the dipoles leading to highest RMS values (Fig G.7). It appears that this subset of dipoles lead to a similar behavior as the entire population. Thus, localization and orientation errors found in this study do not stem from low RMS dipoles only, but most likely from the dipole position and tangentiality in the head, and from geometrical approximation.

5 References

* Algorri, M.-E. and Schmitt, F. Mesh simplification., 1996. In Computer Graphics Forum, volume 15, pages 77–86. Eurographics’96.

* Balish, M., Sato, S., Connaughton, P., and Kufta, C. Localization of implanted dipoles by magnetoencephalography. Neurology, 1991, 41:1072–1076.

* Barth, D., Sutherling, W., Broffman, J., and Beatty, J. magnetic localization of a dipolar current source implanted in a sphere and a human cranium. Electroencephalography and clinical Neurophysiology, 1986, 63:260–273.

* Cohen, D. and Cuffin, B. N. EEG versus MEG localization accuracy: theory and experiment. Brain Topogr, 1991, 4:95–103.

* Cohen, D., Cuffin, B. N., Yunokuchi, K., Maniewski, R., Purcell, C., Cosgrove, G. R., Ives, J., Kennedy, J. G., and Schomer, D. L. MEG versus EEG localization test using implanted sources in the human brain. Ann Neurol, 1990, 28:811–817.

* Cuffin, B. N. Eccentric spheres models of the head. IEEE Trans Biomed Eng, 1991, 38:871–878.

* De Munck, J. C. A linear discretization of the volume conductor boundary integral equation using analytically integrated elements. IEEE Trans Biomed Eng, 1992, 39:986–990.

* Ferguson, A. S., Zhang, X., and Stroink, G. A complete linear discretization for calculating the magnetic field using the boundary element method. IEEE Trans Biomed Eng, 1994, 41:455–460.

* Fuchs, M., Drenckhahn, R., Wischmann, H.-A., and Wagner, M. An improved Boundary Element Method for Realistic Volume-Conductor Modeling. IEEE Trans Biomed Eng, 1998, 45:980–997.
* Hämäläinen, M. Interpretation of neuromagnetic measurements: modeling and statistical considerations. PhD thesis, Epsoo, Hollande, 1987.

* Hämäläinen, M. and Sarvas, J. Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data. IEEE Trans Biomed Eng, 1989, 36:165-171.

* Haueisen, J., Ramon, C., Czapski, P., and Eiselt, M. On the influence of volume currents and extended sources on neuromagnetic fields: a simulation study. Ann Biomed Eng, 1995, 23:728-739.

* Homma, S., Musha, T., Nakajima, Y., Okamoto, Y., Blom, S., Flink, R., and Hagbarth, K. E. Conductivity ratios of the scalp-skull-brain head model in estimating equivalent dipole sources in human brain. Neurosci Res, 1995, 22:51-55.

* Leahy, R., Mosher, J., Spencer, M., Huang, M., and Lewine, J. A study of dipole localization accuracy for MEG and EEG using a human skull phantom. Electroenceph clin Neurophysiol, 1998, 107:159-173.

* Mäkelä, J., Hämäläinen, M., Hari, R., and McEvoy, L. Whole-head mapping of middle-latency auditory evoked magnetic fields. Electroenceph clin Neurophysiol, 1994, 92:414-421.

* Meijs, J. and Peters, M. The EEG and MEG, using a model of eccentric spheres to describe the head. IEEE Trans Biomed Eng, 1987a, 34:913-920.

* Meijs, J., Bosh, F., Peters, M., and Lopes da Silva, F. On the magnetic field distribution generated by a dipolar current source situated in a realistically shaped compartment model of the head. Electroenceph clin Neurophysiol, 1987b, 66:222-226.

* Menninghaus, E. and Lütkenhöner, B. How silent are deep and radial sources in neuromagnetic measurements ?, 1994a. In Proceedings of the Biomag Conference, Vienna.

* Menninghaus, E., Lütkenhöner, B., and Gonzalez, S. L. Localization of a dipolar source in a skull phantom: realistic versus spherical model. IEEE Trans Biomed Eng, 1994b, 41:986-989.

* Merlet, I., Garcia-Larrea, L., Ryvlin, P., Isnard, J., Sindou, M., and Mauguière, F. Topographical reliability of mesio-temporal sources of interictal spikes in temporal lobe epilepsy. Electroenceph clin Neurophysiol, 1998, 107:206-212.

* Rose, D. F., Sato, S., Ducla-Soares, E., and Kufta, C. V. Magnetoencephalographic localization of subdural dipoles in a patient with temporal lobe epilepsy. Epilepsia, 1991, 32:635-641.

* Rush, S. and Driscoll, D. A. Current distribution in the brain from surface electrodes. Anesth Analg, 1968, 47:717-723.

* Scherg, M. Fundamentals of dipole source potential analysis., 1990. In F. Grandori, M. Hoke, G. R., editor, Auditory Evoked Magnetic Fields and Electric Potentials. Adv. Audiol. Basel: Karger, volume 6, pages 40-69.
* Tomita, S., Kajihara, S., Kondo, Y., Yoshida, Y., Shibata, K., and Kado, H. Influence of Head Model in Biomagnetic Source Localization. Brain Topogr., 1996, 8:337-340.

* Verkindt, C., Bertrand, O., Perrin, F., Echallier, J.-F., and Pernier, J. Tonotopic organization of the human auditory cortex: N100 topography and multiple dipole model analysis. Electroenceph clin Neurophysiol, 1995, 96:143-156.

* Yvert, B., Bertrand, O., Echallier, J.-F., and Pernier, J. Improved forward EEG calculations using local mesh refinement of realistic head geometries. Electroencephalogr Clin Neurophysiol, 1995, 95:381-392.

* Yvert, B., Bertrand, O., Echallier, J.-F., and Pernier, J. A systematic evaluation of the spherical model accuracy in EEG dipole localization. Electroenceph clin Neurophysiol, 1997, 102:452-459.

* Yvert, B., Bertrand, O., Pernier, J., and Ilmoniemi, R. Human cortical responses evoked by dichotically presented tones of different frequencies. Neuroreport, 1998, 9:1115-1119.

* Zanow, F. Realistically shaped models of the head and their application to EEG and MEG. PhD thesis, university of Twente, Holland, 1997. pp. 61-75.
<figure><img src="image_92.png" /></figure>
# Annexe H

## Méthode de Choleski

Lorsque j’ai débuté ce travail de thèse, la résolution de la partie linéaire du problème inverse était réalisée en utilisant la décomposition de Choleski (voir section 1.2.4.2). Cette méthode permet de trouver des solutions dans tous les cas, y compris lorsqu’il y a des valeurs propres nulles. Toutefois je souhaitais que la solution proposée dans ce cas soit orthogonale au noyau, pour limiter par exemple une composante radiale en MEG sphérique. Une décomposition en valeur singulière aurait permis de trouver cette solution mais serait plus coûteuse en temps de calcul. C’est pourquoi j’ai implémenté, uniquement dans le cas où apparaissait une valeur propre nulle, une petite variante de l’algorithme de Choleski que j’ai décrite ici.

### 1 Cas général

On veut résoudre le problème

$$SX = Y$$

connaissant $S$ et $Y$. $S$ est une matrice symétrique définie positive. L’algorithme de Choleski consiste à trouver une matrice triangulaire inférieure $L$, telle que $S = LL^T$.

La résolution de notre problème peut alors se faire très simplement en résolvant dans un premier temps :

$$LV = Y$$

et dans un deuxième temps :

$$L^T X = V.$$

<footnote n="0">⁰Cette annexe est référencée à la section 1.2.4.2, page 27</footnote>
Ces deux systèmes sont triangulaires, et donc très faciles à résoudre, le premier en descendant les équations :

$$v_1 = \frac{y_1}{L_{1,1}}$$

et, connaissant $v_k \forall k < i$

$$v_i = \frac{y_i}{L_{i,i}} - \sum_{k=1}^{i-1} \frac{L_{i,k}}{L_{i,i}} v_k$$

et le second en remontant les équations.

Il reste à trouver la matrice L. Pour cela, il suffit d'écrire les équations de $S = LL^T$ pour la partie triangulaire supérieure de S :

$$S_{1,1} = L_{1,1}^2 \qquad (H.1)$$

$$S_{1,i} = L_{1,1} * L_{i,1} \qquad (H.2)$$

$\vdots$

$$S_{i,i} = \sum_{k=1}^{i} L_{i,k}^2 \qquad (H.3)$$

$$S_{i,j} = \sum_{k=1}^{i} L_{i,k} * L_{j,k} \quad j > i \qquad (H.4)$$

On peut donc calculer $L_{i,j}$ en descendant les équations :

$$L_{1,1} = \sqrt{S_{1,1}}$$

$$L_{i,1} = \frac{S_{1,i}}{L_{1,1}}$$

et, connaissant $L_{j,k} \forall k < i, \forall j \ge k$

$$L_{i,i} = \sqrt{S_{i,i} - \sum_{k=1}^{i-1} L_{i,k}^2}$$

$$L_{j,i} = \frac{S_{j,i} - \sum_{k=1}^{i-1} L_{i,k} * L_{j,k}}{L_{i,i}} \quad j > i$$

Ceci permet de calculer $L$ à condition que ce qu'il y a sous les racines carrées soit bien strictement positif. Cette condition nous sera donnée par la positivité de S. En effet,

$$S_{1,1} = \begin{pmatrix} 1 & 0 & \dots & 0 \end{pmatrix} S \begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix} > 0$$
et, en notant $L_{i-1}$ la sous matrice de $L$ constituée de ses $i-1$ premières lignes et colonnes, qui sont déjà connues, et $c$ les $i-1$ premières lignes de la $i$-ème colonne, également déjà connues, on sait que $S$ est de la forme :

$$S = \begin{pmatrix} L_{i-1}L_{i-1}^T & L_{i-1}\mathbf{c}_i & \dots & L_{i-1}\mathbf{c}_j & \dots \\ \mathbf{c}_i^T L_{i-1}^T & S_{i,i} & \dots & S_{i,j} & \vdots \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{c}_j^T L_{i-1}^T & S_{i,j} & \dots & S_{j,j} & \vdots \\ \vdots & \dots & \dots & \dots & \ddots \end{pmatrix}$$

De plus $L_{i-1}$ est inversible puisque $S$ est définie, et donc :

$$S_{i,i} - \sum_{k=1}^{i-1} L_{i,k}^2 = S_{i,i} - \mathbf{c}_i^T \mathbf{c}_i = \left( \mathbf{c}_i^T L_{i-1}^{-1} \quad -1 \quad 0 \quad \dots \quad 0 \right) S \begin{pmatrix} L_{i-1}^{-T} \mathbf{c}_i \\ -1 \\ 0 \\ \vdots \\ 0 \end{pmatrix} > 0$$

## 2 Cas des matrices symétriques positives non définies

Dans le cas des matrices symétriques positives non définies, ce qui se trouve sous les racines carrées reste positif, mais peut être nul ($L_{i,i} = 0$). Dans ce cas, l'équation H.4 se résume à $S_{i,j} = \sum_{k=1}^{i-1} L_{i,k} * L_{j,k} \quad j > i$, et ne donne pas d'information sur $L_{j,i}$. On peut
montrer toutefois, que cette équation est respectée :

$$ \forall \varepsilon, \quad \left( \mathbf{c}_i^T L_{i-1}^{-1} \begin{pmatrix} -1 & 0 & \dots & 0 & \varepsilon & 0 & \dots & 0 \end{pmatrix} \right) S \begin{pmatrix} L_{i-1}^{-T} \mathbf{c}_i \\ -1 \\ 0 \\ \vdots \\ 0 \\ \varepsilon \\ 0 \\ \vdots \\ 0 \end{pmatrix} \geq 0 $$

$$ \forall \varepsilon, \quad (S_{i,i} - \mathbf{c}_i^T \mathbf{c}_i) + \varepsilon(2(\mathbf{c}_i^T \mathbf{c}_j - S_{i,j}) + \varepsilon S_{j,j}) \geq 0 $$

$$ \forall \varepsilon, \quad L_{i,i}^2 + \varepsilon(2(\mathbf{c}_i^T \mathbf{c}_j - S_{i,j}) + \varepsilon S_{j,j}) \geq 0 $$

$$ \forall \varepsilon, \quad \varepsilon(2(\mathbf{c}_i^T \mathbf{c}_j - S_{i,j}) + \varepsilon S_{j,j}) \geq 0 $$

$$ \left\{ \begin{array}{l} \lim_{\varepsilon \to 0^+} (2(\mathbf{c}_i^T \mathbf{c}_j - S_{i,j}) + \varepsilon S_{j,j}) \geq 0 \\ \lim_{\varepsilon \to 0^-} (2(\mathbf{c}_i^T \mathbf{c}_j - S_{i,j}) + \varepsilon S_{j,j}) \leq 0 \end{array} \right. $$

$$ \left\{ \begin{array}{l} \mathbf{c}_i^T \mathbf{c}_j - S_{i,j} \geq 0 \\ \mathbf{c}_i^T \mathbf{c}_j - S_{i,j} \leq 0 \\ \mathbf{c}_i^T \mathbf{c}_j - S_{i,j} = 0 \end{array} \right. $$

$$ S_{i,j} = \sum_{k=1}^{i-1} L_{i,k} * L_{j,k} $$

On peut donc choisir arbitrairement $L_{j,i} = 0$ dans ce cas et trouver ainsi une factorisation de Choleski de notre matrice.

Les systèmes linéaires

$$ SX = Y $$

ont alors soit plusieurs solutions, soit aucune selon que $Y$ appartienne à l'espace image de $S$ ou non. Notre but peut être soit de fournir une de ces solutions, soit de fournir la solution qui est orthogonale au noyau de $S$. Dans tous les cas, on décompose le problème en deux sous-problèmes :

$$ \begin{cases} LV = Y \\ L^T X = V \end{cases} $$
Comme dans le cas des matrices définies, on résout ces problèmes par une descente et une remontée :

$$
\begin{array}{ll}
\text{descente :} &
\left\{
\begin{array}{l}
V_1 = \frac{Y_1}{L_{1,1}} \text{ si } L_{1,1} \neq 0 \\
V_1 = \mathbf{c}_1^T \mathbf{X} = 0 \text{ si } L_{1,1} = 0 \\
\text{et, connaissant } V_k \quad \forall k < i \\
V_i = \frac{Y_i}{L_{i,i}} - \sum_{k=1}^{i-1} \frac{L_{i,k}}{L_{i,i}} V_k \text{ si } L_{i,i} \neq 0 \\
V_i = \mathbf{c}_i^T \mathbf{X} = 0 \text{ si } L_{i,i} = 0
\end{array}
\right. \\
\\
\text{remontée :} &
\left\{
\begin{array}{l}
X_n = \frac{V_n}{L_{n,n}} \text{ si } L_{n,n} \neq 0 \\
X_n = X_n \text{ Peut être fixé arbitrairement si } L_{1,1} = 0 \\
\text{et, connaissant } X_k \quad \forall k > i \\
X_i = \frac{V_i}{L_{i,i}} - \sum_{k=i+1}^{n} \frac{L_{k,i}}{L_{i,i}} X_k \text{ si } L_{i,i} \neq 0 \\
X_i = X_i \text{ Peut être fixé arbitrairement si } L_{i,i} = 0
\end{array}
\right.
\end{array}
$$

De plus, des conditions apparaissent à l'existence d'une solution au système :

$$ \text{si } L_{i,i} = 0, Y_i = \sum_{\substack{k=1 \\ L_{k,k} \neq 0}}^{i-1} L_{i,k} V_k $$

Si ces conditions sont réunies, pour n'importe quelles valeurs des {$X_i, L_{i,i} = 0$}, on obtient une solution. Le choix de Nash [1979] est de les choisir égaux à 0, on obtient ainsi une solution que nous noterons par la suite $\tilde{X}$. Toutefois, on peut également souhaiter trouver les solutions orthogonales au noyau, par exemple dans le cas de la MEG sphérique, lorsque le noyau correspond aux dipôles radiaux.

Dans ce cas, le problème consiste à trouver $X$ tel que :

$$
\left\{
\begin{array}{l}
\forall \mathbf{Z}, LL^T \mathbf{Z} = 0 \Rightarrow \mathbf{X}^T \mathbf{Z} = 0 \\
L^T \mathbf{X} = \mathbf{V}
\end{array}
\right.
$$

or $LL^T \mathbf{Z} = 0 \Rightarrow \mathbf{Z}^T LL^T \mathbf{Z} = 0 \Rightarrow \|L^T \mathbf{Z}\| = 0 \Rightarrow L^T \mathbf{Z} = 0$, donc notre problème équivaut à :

$$
\left\{
\begin{array}{l}
\forall \mathbf{Z}, L^T \mathbf{Z} = 0 \Rightarrow \mathbf{X}^T \mathbf{Z} = 0 \\
L^T \mathbf{X} = \mathbf{V}
\end{array}
\right.
$$

En choisissant $\mathbf{X}$ sous la forme $\mathbf{X} = L\lambda$, le problème devient :

$$
\left\{
\begin{array}{l}
\forall \mathbf{Z}, L^T \mathbf{Z} = 0 \Rightarrow \lambda^T L^T \mathbf{Z} = 0 \\
L^T L \lambda = \mathbf{V}
\end{array}
\right.
$$
et donc plus simplement $L^T L \lambda = V$. L'existence d'au moins une solution à ce problème est garantie par le fait que $\exists \tilde{\mathbf{X}}, \mathbf{V} = L^T \tilde{\mathbf{X}}$. En effet, le problème équivaut à :

$$L^T L \lambda = L^T \tilde{\mathbf{X}}$$

or cette classe de problème a toujours au moins une solution, et ce, quelle que soit la matrice $L$, triangulaire ou non. En effet, la matrice $L^T L$ est symétrique, donc elle est diagonalisable, et donc il existe une matrice de changement de base $P$ telle que $P^T L^T L P$ soit diagonale à diagonale croissante :

$$P^T L^T L P = \begin{pmatrix} 0 & & & & \\ & \ddots & & & \\ & & \mathbf{0} & & \\ & & & d_1 & \\ \mathbf{0} & & & & \ddots \\ & & & & d_p \end{pmatrix}$$

Et comme $P$ est inversible d'inverse $P^T$, notre problème équivaut à :

$$P^T L^T L P P^T \lambda = P^T L^T \tilde{\mathbf{X}}.$$

Posons $\mathbf{U} = P^T \lambda$, et $\mathbf{W} = P^T L^T \tilde{\mathbf{X}}$. Notre problème équivaut à résoudre :

$$\begin{pmatrix} 0 & & & & \\ & \ddots & & & \\ & & \mathbf{0} & & \\ & & & d_1 & \\ \mathbf{0} & & & & \ddots \\ & & & & d_p \end{pmatrix} \begin{pmatrix} \mathbf{U}_1 \\ \vdots \\ \mathbf{U}_{n-p} \\ \mathbf{U}_{n-p+1} \\ \vdots \\ \mathbf{U}_n \end{pmatrix} = \begin{pmatrix} \mathbf{W}_1 \\ \vdots \\ \mathbf{W}_{n-p} \\ \mathbf{W}_{n-p+1} \\ \vdots \\ \mathbf{W}_n \end{pmatrix}$$

Ce qui équivaut à

$$\begin{cases} \forall i \in [1..n-p], & \mathbf{W}_i = 0 \\ \forall i \in [n-p+1..n], & \mathbf{U}_i = \frac{1}{d_{i+p-n}} \mathbf{W}_i \end{cases}$$

Les premières équations sont les équations de consistance de notre problème. Dans notre cas, elles sont toujours vérifiées. En effet, posons $Z = LP$. La nullité des $n-p$ premiers termes diagonaux implique :

$$\forall i \in [1..n-p], \sum_{k=1}^{n} Z_{k,i}^2 = 0$$
et donc,

$$ \forall i \in [1..n-p], \forall k \in [1..n], Z_{k,i} = 0 $$

Or $\mathbf{W} = Z^T \tilde{\mathbf{X}}$. Donc :

$$ \forall i \in [1..n-p], W_i = \sum_{k=1}^{n} Z_{k,i} \tilde{X}_k = 0 $$

En conclusion, le vecteur $\mathbf{X} = L\lambda$ tel que $\lambda$ est solution de $L^T L\lambda = V$ est la solution de $SX = Y$ orthogonale au noyau de $S$.
<figure><img src="image_93.png" /></figure>
# Annexe I

## Génération d'un bruit réaliste

Pour générer un bruit présentant une corrélation spatiale et intermodale (entre les voies MEG et les électrodes EEG) réaliste, nous avons utilisé un enregistrement d'EEG et de MEG réalisé sur un patient épileptique alors qu'on ne lui demandait aucune activité particulière, et ne présentant pas de pointes épileptiques. Cet enregistrement a été réalisé au centre MEG de l'hôpital de la Pitié Salpêtrière à PARIS sur une MEG CTF à 151 canaux de type gradiomètre du troisième ordre, et simultanément sur 32 voies d'EEG, et nous a été fourni gracieusement par Line GARNERO. La fréquence d'échantillonage était de 625 Hz, et lors de l'enregistrement un filtrage passe-bas de type butterworth avec une fréquence de coupure à 200 Hz a été appliqué, ainsi que des filtrages coupe-bande de type butterworth (notch) à 50 Hz avec une largeur de bande de 5 Hz (appliqué deux fois), à 100 Hz avec une largeur de bande de 5 Hz, et à 150 Hz avec une largeur de bande de 10 Hz.

### 1 Traitement d'un enregistrement couplé spontané pour générer un bruit réaliste

#### 1.1 Recalage des électrodes et des canaux MEG sur un modèle réaliste

L'enregistrement utilisé pour générer le bruit réaliste provenait d'un sujet différent de celui que nous avions utilisé pour construire notre modèle réaliste. Nous avons donc du recaler au mieux les électrodes et les capteurs MEG sur notre modèle réaliste.

Les positions des voies EEG (voir figure I.1) ont été utilisées pour déterminer la position de la tête du sujet afin de mettre les coordonnées des voies (qui étaient en repère constructeur) dans le repère lié à la tête utilisé dans notre unité. En effet, le modèle réaliste que nous avons utilisé a été obtenu chez un sujet différent et est décrit dans ce repère. Ce repère (qui sera

<footnote n="0">⁰Cette annexe est référencée à la section 4.4.1, page 99</footnote>
<figure><img src="image_94.png" /><figcaption>FIG. I.1: Positions des voies EEG.</figcaption></figure>

<figure><img src="image_95.png" /><figcaption>FIG. I.2: Détermination du repère 10-20 à partir des positions des électrodes.</figcaption></figure>

nommé ensuite repère tête de l'unité 280) est défini comme suit : l'axe x passe par les tragus, l'axe y est la perpendiculaire à l'axe x qui passe par le nasion, et l'axe z est perpendiculaire aux axes x et y.

Il faut donc déterminer les positions approximatives des tragus et du nasion dans le repère du constructeur pour pouvoir trouver la matrice de passage du repère constructeur au repère tête de l'unité 280. Pour cela, nous avons utilisé leurs coordonnées habituelles dans le système 10-20.
<figure><img src="image_96.png" /><figcaption>FIG. I.3: Passage du repère 10-20 au repère tête de l'unité 280.</figcaption></figure>

Les milieux des électrodes (FP1, FP2), (F7, F8), (T7, T8), et (TP7, TP8) ainsi que les électrodes FPZ et OZ devraient théoriquement, dans le système 10-20, se situer sur l'axe y du système 10-20. Ils nous permettent de définir un axe y moyen pour le repère 10-20. Le centre du système 10-20 est déterminé par l'intersection de cet axe avec le plan moyen passant par T3 (ici équivalent à T7), C3, Cz, C4 et T4 (ici équivalent à T8). La direction yz est alors donnée par la droite partant du centre et passant par le barycentre de C3,C4 et Cz. L'axe x du système 10-20 est orthogonal aux axes y et à la direction yz, et l'axe z peut alors être trouvé comme la perpendiculaire aux axes x et y. Nous pouvons maintenant passer du repère 10-20 au repère constructeur. Nous utilisons les coordonnées habituelles des tragus et du nasion dans le système 10-20 : en coordonnées sphériques $(\theta(°), \phi(°))$ : Nasion=(110, 90.0), tragus gauche=(120, 180.0), tragus droit=(120, 0.0) pour déterminer les coordonnées de ces points dans le repère constructeur. Ceux-ci nous permettent alors de déterminer le changement de repère à appliquer à toutes les coordonnées des voies (MEG et EEG) pour passer dans le repère tête de l'unité 280. Les électrodes ont ensuite été projetées par projection centrale sur notre modèle de tête (décrit à la section 4.5) pour se trouver sur le scalp de notre modèle réaliste (voir figure I.4).

## 1.2 Traitement des signaux

Trois des 32 voies d'EEG ont dû être exclues soit parce que leur coordonnées spatiales étaient inconnues soit parce qu'elles présentaient un niveau de bruit trop élevé. Il restait donc 151 voies MEG et 29 voies EEG. Les fichiers contenant les enregistrements ont été convertis au format "eeg" de l'unité 280 afin de pouvoir être traités. Comme des harmoniques de 50 Hz persistaient dans le spectre à des fréquences supérieures à 200 Hz, un filtre passe bas avec une
<figure><img src="image_97.png" /><figcaption>FIG. I.4: Positions des voies MEG (en rouge) et EEG (en bleu) recalées sur notre modèle.</figcaption></figure>

fréquence de coupure à 150 Hz a été appliqué (voir figure I.5). Ce filtrage, comme le filtrage analogique d'origine, crée une certaine corrélation temporelle.

<figure><img src="image_98.png" /><figcaption>FIG. I.5: Spectres d'une voie MEG (à gauche) et d'une voie d'EEG (à droite).</figcaption></figure>

Les enregistrements spontanés MEG et EEG ont ensuite été moyennés afin de reproduire de manière réaliste les conditions d'une réponse évoquée. Ce moyennage peut avoir pour conséquence de réduire la corrélation spatiale. Afin de se mettre dans des situations de bruit ayant des corrélations différentes, deux façons de moyenner ont été utilisées : dans un cas, (bruit1) la durée retenue pour chaque intervalle de temps était de 200 échantillons, et les
échantillons de départ étaient pris tous les 100 + rand40 échantillons, où rand40 est un entier aléatoire uniformément distribué entre 0 et 40. On a ainsi obtenu 200 réalisations de bruit résultant de la moyenne de 560 essais élémentaires. Dans le second cas, (bruit2) les échantillons de départ étaient pris tous les 50 + rand20 échantillons, et 100 réalisations de bruit résultant de la moyenne de 1184 essais élémentaires ont été obtenus. La fréquence d'échantillonage était de 625 échantillons par seconde. Les enregistrements contenaient 40 séquences de 1875 échantillons. Les intervalles contenant des ruptures d'acquisition ont été exclus dans les moyennes. Un centrage des valeurs moyennées a été appliqué à chaque voie. On a ainsi obtenu pour chaque voie 200 valeurs d'un bruit de moyenne nulle, et dont les écarts-type sont cartographiés figure I.6. Une analyse en composantes principale des réalisations de bruit ainsi obtenus donne comme première composante principale celle représentée figure I.7. Quelques topographies de bruit sont présentées figure I.8. Enfin, l'intercorrélation spatiale entre les voies est représentée graphiquement figure I.9. On constate que cette intercorrélation peut atteindre des niveaux assez important (carrés rouges) y compris entre les différentes modalités (rectangles en bas à droite et en haut à gauche).

<figure><img src="image_99.png" /><figcaption>FIG. I.6: Cartes des écarts types par voie des réalisations de bruit, pour l'EEG (à gauche) et la MEG (à droite), bruit1 (en haut) et bruit2 (en bas).</figcaption></figure>

À partir des réalisations de bruit ainsi obtenues pour 29 électrodes, une interpolation spline sphérique a été utilisée pour obtenir des réalisations de bruit pour des configurations à 64, 151 et 256 électrodes représentées figure I.10.
<figure><img src="image_100.png" /><figcaption>FIG. I.7: Cartes de la première composante principale, pour l'EEG (à gauche) et la MEG (à droite), bruit1 (en haut) et bruit2 (en bas). La première composante principale explique 72.5 % de la variance pour bruit1 en EEG, 39.5 % pour bruit2 en EEG, 39.7 % de la variance pour bruit1 en MEG, et 60.7 % pour bruit2 en MEG</figcaption></figure>
<figure><img src="image_101.png" /><figcaption>FIG. I.8: Topographies de quatre réalisations de bruit. De haut en bas : eeg, bruit1 ; eeg, bruit2 ; meg, bruit1 ; meg, bruit2.</figcaption></figure>

<figure><img src="image_102.png" /><figcaption>FIG. I.9: Représentation graphique de l'intercorrélation normalisée entre les voies. Les 29 premières voies sont des voies d'EEG (carré en bas à gauche), les 151 dernières voies sont des voies de MEG (carré en haut à droite). Pour une abscisse i et une ordonnée j, la couleur représente la valeur absolue de l'intercorrélation normalisée $|\sigma_{ij}|$ entre la voie i et la voie j. La diagonale ascendante correspond à l'autocorrélation normalisée, qui vaut 1.</figcaption></figure>
<figure><img src="image_103.png" /><figcaption>FIG. I.10: Positions des électrodes dans les configurations à 29, 64, 151 ou 256 électrodes.</figcaption></figure>
# Annexe J

## Validation d'une méthode rapide de calcul du problème direct

Cette annexe présente une méthode d'accélération de la modélisation réaliste en MEG et en EEG utilisant une interpolation des champs engendrés. Dans un premier temps, cette annexe présente la méthode et son évaluation dans un modèle sphérique telle qu'elle a été proposée dans un article en cours de soumission. Ce travail a été réalisé pour la majeure partie par Blaise Yvert et sera présenté plus en détail dans un article en cours de publication [Yvert et al., 2001b]. Dans un deuxième temps, une pré-validation pour les géométries réalistes est présentée. Cette pré-validation a permis l'utilisation de cette méthode d'accélération pour évaluer l'apport du couplage au chapitre 4.

### 1 Présentation et évaluation de la méthode dans le cas d'une géométrie sphérique

#### 1.1 Introduction

L'utilisation de modèles réalistes pour résoudre les problèmes direct et inverses en EEG/MEG est limitée en pratique par le temps de calcul important des méthodes numériques utilisées. Un compromis doit souvent être recherché entre la précision souhaitée pour le calcul et un temps de calcul acceptable. Nous proposons ici une méthode pour accélérer ces calculs lorsque de nombreux problèmes directs doivent être résolus avec un même modèle réaliste. Cette méthode est compatible avec toutes les méthodes numériques (intégrales de surface, éléments finis, différences finies). Elle peut également être étendue à l'ECG/EMG. Afin

<footnote n="0">⁰Cette annexe est référencée à la section 4.6.1, page 100</footnote>
d'évaluer la précision de cette méthode, un modèle à géométrie sphérique a été utilisé car on connaît la solution analytique du problème direct dans ce modèle.

## 1.2 Méthodes

Pour une position donnée $\mathbf{r}$ dans le cerveau, appelons $\mathbf{L}_{x,0}(\mathbf{r})$, $\mathbf{L}_{y,0}(\mathbf{r})$, et $\mathbf{L}_{z,0}(\mathbf{r})$ les trois champs produits par les dipôles unitaires à la position $\mathbf{r}$ ayant une orientation selon les directions positives $x$, $y$, et $z$ du repère cartésien. Ils forment une base de l'espace engendré par les solutions au problème direct pour tous les dipôles situés à cette position. En effet, tout dipôle situé en $\mathbf{r}$ avec un moment dipolaire

$$ \mathbf{Q} = \begin{bmatrix} Q_x \\ Q_y \\ Q_z \end{bmatrix} $$

produit un champ $\mathbf{L}(\mathbf{r}) = \{V(\mathbf{r}), \mathbf{B}(\mathbf{r})\}$ tel que :

$$ \mathbf{L}(\mathbf{r}) = Q_x \mathbf{L}_{x,0}(\mathbf{r}) + Q_y \mathbf{L}_{y,0}(\mathbf{r}) + Q_z \mathbf{L}_{z,0}(\mathbf{r}) \quad (J.1) $$

Où $V(\mathbf{r})$ et $\mathbf{B}(\mathbf{r})$ sont des vecteurs dont les coefficients sont respectivement le potentiel et le champ magnétique produit à chaque voie par le dipôle de moment $\mathbf{Q}$ en $\mathbf{r}$.

### 1.2.1 approche classique

Lorsqu'on utilise des modèles réalistes et la méthode des intégrales de surface, pour calculer $V(\mathbf{r})$ et $\mathbf{B}(\mathbf{r})$, il faut résoudre les systèmes linéaires suivants :

$$ \mathbf{V} = \mathbf{A}_{eeg}^{-1} \mathbf{V}_o $$

$$ \mathbf{B} = \mathbf{B}_o + \mathbf{A}_{meg} \mathbf{V} $$

où :

* $V$ est le vecteur dont les coefficients sont le potentiel à chacun des $N$ nœuds du maillage,

* $\mathbf{A}_{eeg}$ est une matrice $N \times N$,

* $\mathbf{V}_o$ est un vecteur dépendant de la source et de la géométrie du modèle,

* $\mathbf{B}$ est un vecteur dont les coefficients sont la composante du champ magnétique selon la normale de chacune des $N_c$ spires du capteur.

* $\mathbf{B}_o$ est un vecteur dont les coefficients sont les composantes du champ magnétique aux spires engendrées par les courants primaires uniquement. Il peut être calculé analytiquement.
• $A_{\text{meg}}$ est une matrice $N_c * N$ ($A_{\text{meg}}V$ est la contribution des courants volumiques).

Les matrices $A_{\text{eeg}}$ et $A_{\text{meg}}$ peuvent être calculées une fois pour toutes car elles ne dépendent que de la géométrie et des conductivités des différents compartiments du modèle de la tête, et non des sources. On calcule généralement l'inverse ou la décomposition LU de $A_{\text{eeg}}$ une fois pour toutes. Chaque résolution du premier système linéaire demande encore $N^2$ multiplications et $N^2$ additions, ce qui rend les temps de calcul importants et directement liés à la taille du maillage.

### 1.2.2 approche proposée

**Calculs préliminaires** Tout comme dans l'approche classique, il faut calculer les matrices $A_{\text{eeg}}$, $A_{\text{meg}}$ et l'inverse ou la décomposition LU de $A_{\text{eeg}}$. L'approche proposée consiste à faire un calcul préliminaire supplémentaire. Il suffit de faire ce calcul une seule fois pour un sujet et une configuration de capteurs donnés. Ce calcul supplémentaire consiste en :

1. La construction d'une grille discrète couvrant le volume du cerveau.

2. Le calcul des champs engendrés $L_{x,0}(\mathbf{r})$, $L_{y,0}(\mathbf{r})$, et $L_{z,0}(\mathbf{r})$ en chaque point $\mathbf{r}$ de la grille.

**Résolution approchée du problème direct pour une source quelconque.** Une fois cette étape préliminaire accomplie, chaque problème direct est résolu de manière approchée à partir des champs engendrés pré-calculés. Pour toute source dipolaire $A$ située dans le volume cérébral à la position $\mathbf{R}_a$ avec un moment dipolaire $\mathbf{Q}_a$, le champ $L_a = \{V_a, B_a\}$ engendré par cette source est estimé comme suit : dans un premier temps, les trois champs engendrés $L_{x,0}(\mathbf{R}_a)$, $L_{y,0}(\mathbf{R}_a)$, et $L_{z,0}(\mathbf{R}_a)$ sont interpolés à partir des champs engendrés correspondants aux points de la grille dans le voisinage de $\mathbf{R}_a$. Dans un second temps, le champ $L_a$ correspondant à l'orientation et au moment dipolaire $\mathbf{Q}_a$ est calculé selon la formule de l'équation J.1. Notons que si $A$ est une source périphérique en dehors de la grille, le calcul se fait par extrapolation. Trois techniques d'interpolation ont été utilisées pour cette évaluation : une interpolation trilinéaire, une interpolation de Bézier (qui correspond à une extension de l'interpolation trilinéaire avec des polynômes de degré supérieur), et une interpolation par des splines 3D [Duchon, 1976; Babiloni et al., 1996]. Nous avons utilisé des splines 3D d'ordre 3 avec les 80 points de la grille les plus proches de la source.

### 1.2.3 Procédure d'évaluation

Ces simulations utilisent des modèles de la tête à géométrie sphérique afin de pouvoir
comparer les résultats trouvés avec la solution analytique, ainsi que des grilles isotropiques
<figure><img src="image_104.png" /><figcaption>FIG. J.1: A : Montage d'électrodes et de voies MEG utilisées dans nos simulations. B : Grille utilisée pour le précalcul des champs engendrés et grille fine utilisée pour l'évaluation. C : Positions des sources pour le problème inverse.</figcaption></figure>

avec 6 pas de grille différents : 6, 8, 10, 12, 15 et 20 mm, ce qui donne respectivement 9315, 3884, 2007, 1189, 619 et 251 points dans la grille. Dans un premier temps, pour rendre compte de l'erreur intrinsèque due à l'interpolation, la méthode a été testée avec des champs engendrés calculés analytiquement en chaque point de la grille. Dans un deuxième temps, les champs engendrés ont été calculé avec la méthode des intégrales de surface et l'approche LINT [De Munck, 1992; Ferguson et al., 1994]. La méthode a été évaluée pour les problèmes direct et inverse en MEG et en EEG. La figure J.1-A montre les 63 électrodes EEG, le montage MEG à 143 canaux, et le modèle à 1500 nœuds par couche. Le problème direct a été évalué en chaque point d'une grille plus fine dont le pas est de 4.3 mm. En EEG, trois orientations du dipôle unitaire sont considérées, en MEG, deux orientations tangentielles selon les directions $\theta$ et $\phi$ du repère sphérique local. Les solutions trouvées par interpolation sont comparées aux solutions calculées analytiquement selon la formule :

$$RMS = \sqrt{\frac{\sum_{ori=1}^{3} \sum_{v=1}^{Nbvoies} (l_{v,ori}^{ana} - l_{v,ori}^{interp})^2}{\sum_{ori=1}^{3} \sum_{v=1}^{Nbvoies} (l_{v,ori}^{ana})^2}}$$

où $l_{v,ori}$ est le potentiel ou le champ magnétique à la voie v pour l'orientation ori. La figure J.1-B donne un exemple de grille utilisée pour les précalculs de champs engendrés (ronds noirs) et de grille fine où les champs engendrés sont calculés par interpolation (petites croix noires). Pour évaluer le problème inverse, nous avons utilisé 261 positions de dipôles sur 9 axes d'un quadrant supérieur de la sphère (voir figure J.1-C) avec 3 orientations selon les axes x, y et z pour l'EEG (soit 783 dipôles au total) et 2 orientations tangentielles pour la MEG, soit 522 dipôles. Pour chaque dipôle, le problème direct a été résolu analytiquement pour générer les données simulées, et le problème inverse a été résolu sur ces données simulées en utilisant soit la méthode des intégrales de surface, soit notre méthode par interpolation à partir d'une grille pour laquelle les champs engendrés sont pré-calculés soit
analytiquement, soit par la méthode des intégrales de surface. En MEG, comme le modèle réaliste utilisé est un maillage de sphère, on a contraint l'orientation du dipôle à rester tangentielle. Les erreurs d'orientation et de localisation entre les dipôles trouvés par résolution du problème inverse et les dipôles d'origine ont été moyennées entre les axes et les différentes orientations de façon à obtenir une erreur moyenne pour chaque profondeur du dipôle.

### 1.2.4 Effet de troncature de la grille

Pour les positions de sources les plus périphériques, l'erreur numérique peut devenir importante lors de la résolution du problème direct, en particulier en EEG lorsqu'aucun affinage du maillage n'est utilisé [Yvert et al., 1995; 1996]. Cependant, l'utilisation de modèles affinés localement n'est pas facilement applicable en pratique car le modèle devrait être modifié au fur et à mesure des variations de la position de la source lors de la résolution du problème inverse, ce qui suppose à chaque fois un nouveau calcul des matrices $A_{eeg}$ et $A_{meg}$ et une nouvelle inversion de $A_{eeg}$. Intuitivement, afin d'éviter d'introduire une erreur numérique trop importante, les points les plus périphériques de la grille devraient être écartés. Nous avons testé si des grilles tronquées permettraient d'améliorer la précision de l'interpolation en écartant les points de la grille dont la profondeur est inférieure à la longueur moyenne des triangles (8 mm). Rappelons que si une source périphérique se trouve en dehors de la grille, le calcul se fait par extrapolation.

## 1.3 Résultats

### 1.3.1 Temps de calcul

Les calculs préliminaires correspondant au calcul des matrices $A_{eeg}$ et $A_{meg}$ et à l'inversion ou la décomposition LU de $A_{eeg}$ sont communs à l'approche classique et à l'approche proposée. Dans cette simulation, il a fallu 61 minutes à un Pentium-II cadencé à 333 MHz avec 256 Mo de RAM pour faire ces calculs (en utilisant la décomposition LU de $A_{eeg}$). L'approche proposée entraîne une étape de calculs préliminaires supplémentaires dont la durée est proportionnelle au nombre de points dans la grille. La table J.1 donne ces temps de calcul pour chaque pas de grille. Elle montre également qu'un modèle réaliste compte environ 20% de points en moins qu'un modèle sphérique de rayon 78.3 mm pour un même pas de grille, ce qui réduit d'autant le temps de calcul.

La Figure J.2 montre les temps nécessaires pour calculer les champs engendrés par chaque point de la grille fine de pas 4.3 mm avec les différentes approches. Comme on pourrait s'y attendre, la méthode la plus rapide est l'interpolation trilinéaire qui permet de diviser le temps de calcul par environ 20000 en EEG et 13700 en MEG par rapport à l'approche classique.
<table><caption>TAB. J.1: Temps de calcul des champs engendrés (étape préliminaire de la méthode proposée) pour chaque grille utilisée dans cette étude. La dernière colonne montre le nombre de points de la grille qui serait obtenu dans un modèle réaliste. On constate qu'une grille réaliste contient environ 20% de points en moins qu'une sphère de rayon 78.3 mm.</caption><thead><tr><th>type de grille</th><th>pas (mm)</th><th>nombre de points de la grille</th><th>temps de calcul en EEG</th><th>temps de calcul en MEG</th><th>nombre de points dans un modèle réaliste</th></tr></thead><tbody><tr><td>non tronquées</td><td>4.3</td><td>25341</td><td>28h08'</td><td>34h27'</td><td>20183</td></tr><tr><td>tronquées</td><td>6</td><td>9315</td><td>10h07'</td><td>12h16'</td><td>7444</td></tr><tr><td></td><td>8</td><td>3887</td><td>4h18'</td><td>5h05'</td><td>3144</td></tr><tr><td></td><td>10</td><td>2007</td><td>2h14'</td><td>2h34'</td><td>1608</td></tr><tr><td></td><td>12</td><td>1189</td><td>1h18'</td><td>1h31'</td><td>944</td></tr><tr><td></td><td>15</td><td>619</td><td>40'</td><td>47'</td><td>475</td></tr><tr><td></td><td>20</td><td>251</td><td>16'</td><td>19'</td><td>199</td></tr><tr><td>tronquées</td><td>6</td><td>6667</td><td>7h21'</td><td>8h40'</td><td>5156</td></tr><tr><td></td><td>8</td><td>2801</td><td>3h06'</td><td>3h40'</td><td>2161</td></tr><tr><td></td><td>10</td><td>1365</td><td>1h31'</td><td>1h47'</td><td>1110</td></tr><tr><td></td><td>12</td><td>847</td><td>54'</td><td>1h06'</td><td>635</td></tr><tr><td></td><td>15</td><td>437</td><td>28'</td><td>36'</td><td>330</td></tr><tr><td></td><td>20</td><td>179</td><td>11'</td><td>15'</td><td>132</td></tr></tbody></table>

<figure><img src="image_105.png" /><figcaption>FIG. J.2: Temps de calcul pour les différentes approches.</figcaption></figure>
L'interpolation de Bézier permet de gagner un facteur environ 1200 en EEG et 750 en MEG, et les splines 3D un facteur d'environ 40 en EEG et 23 en MEG. Alors que le temps de calcul avec la méthode des intégrales de surfaces dépend davantage de la taille du maillage que du nombre de voies, le temps de calcul pour l'approche proposée est proportionnel au nombre de voies. Ainsi, l'approche proposée est d'autant plus efficace en terme d'accélération des calculs que le nombre de voies est faible. Mais même pour un nombre de voies important comme dans notre simulation, l'accélération est remarquable.

### 1.3.2 Précision intrinsèque de la méthode (précalculs analytiques)

Pour évaluer la précision intrinsèque de la méthode, nous avons effectué les précalculs en chaque point de la grille de manière analytique dans le cas sphérique. Ainsi, les erreurs ne peuvent pas être attribuées à l'approximation numérique. Les erreurs RMS pour l'EEG et la MEG lors de la résolution du problème direct sont représentées figure J.3-A. Les erreurs de localisation et d'orientation lors de la résolution du problème inverse sont représentées figure J.3-B. Ces courbes montrent qu'une précision intrinsèque inférieure à 0.5 mm et 0.2° peut être obtenue pour toutes les profondeurs de dipôles avec une grille suffisamment fine (pas de 6 à 8 mm) et une interpolation Bézier ou Spline 3D. Lorsque les dipôles sont profonds, des grilles plus grossières suffisent à assurer la même précision. L'interpolation trilinéaire, bien que plus rapide que les deux autres, donne des résultats moins précis, en particulier en MEG, et en particulier pour les dipôles profonds.

### 1.3.3 Précision de la méthode (précalculs par intégrales de surface)

Cette fois-ci, pour des grilles non tronquées, les précalculs sont effectués en chaque point de la grille avec la méthode des intégrales de surface et notre maillage de sphère. Les problèmes directs et inverses sont résolus avec la méthode proposée, et les erreurs sont comparées à celles obtenues avec la méthode des intégrales de surface classique (Figure J.4). On y constate que les erreurs sont très comparables à l'approche classique. Les erreurs dûes à l'interpolation sont négligeables par rapport aux erreurs numériques. Pour les dipôles dont la profondeur est supérieure à 20 mm, des grilles de pas 10-12 mm suffisent à obtenir une bonne précision.

### 1.3.4 Effet de troncature

La figure J.5 montre les erreurs obtenues pour des grilles tronquées. Elle montre une amélioration de l'erreur RMS grâce à la troncature pour les dipôles périphériques en EEG. L'erreur d'orientation est également améliorée en EEG.
A)
Forward
Problem

<figure><img src="image_106.png" /></figure>

B) Inverse Problem

<figure><img src="image_107.png" /><figcaption>FIG. J.3: Erreurs intrinsèques à la méthode en fonction de la coordonnée radiale du dipôle (en abscisse). A : Problème direct. B : Problème inverse, erreurs de localisation (à gauche), et d'orientation (à droite).</figcaption></figure>
A) Forward Problem

<figure><img src="image_108.png" /></figure>

B) Inverse Problem

<figure><img src="image_109.png" /><figcaption>FIG. J.4: Grilles non tronquées. Erreurs obtenues en fonction de la coordonnée radiale du dipôle (en abscisse). A : Problème direct. B : Problème inverse, erreurs de localisation (à gauche), et d'orientation (à droite). BEM = méthode des intégrales de surface.</figcaption></figure>
A) Forward Problem

<figure><img src="image_110.png" /></figure>

B) Inverse Problem

<figure><img src="image_111.png" /></figure>

<figure><img src="image_112.png" /><figcaption>FIG. J.5: Grilles tronquées. Erreurs obtenues en fonction de la coordonnée radiale du dipôle (en abscisse). A : Problème direct. B : Problème inverse, erreurs de localisation (à gauche), et d'orientation (à droite). BEM = méthode des intégrales de surface.</figcaption></figure>
## 2 Validation préliminaire dans le cas d'une géométrie réaliste en MEG

Pour valider la possibilité d'utiliser les grilles en MEG, modalité pour laquelle la géométrie sphérique représente un cas vraiment particulier, j'ai effectué une évaluation par simulation avec les 2000 dipôles qui avaient servi pour le chapitre 3.

Le but n'était pas de faire une étude approfondie des différences entre l'utilisation de la grille et des intégrales de surface classiques, mais de vérifier que les grilles semblaient bien faire leur office et pourraient donc être utilisées. En particulier, leur application pour calculer des intervalles de confiance avec un modèle réaliste a été étudiée.

Comme pour l'article d'évaluation de la précision du problème inverse en MEG, les données étaient des données simulées de champ magnétique calculées avec la méthode des intégrales de surface pour un maillage très fin (3L1500). Le problème inverse a alors été résolu en utilisant soit les grilles, soit la méthode habituelle, avec un modèle réaliste plus grossier (3L900).

### 2.1 Méthodes

Nous avons utilisé une grille de pas 10 mm tronquée à 2 mm du maillage. Le temps de calcul pour les champs résultants a été de 34 minutes sur un pentiumIII à 450 MHz. Une interpolation spline 3D a été utilisée pour estimer le champ magnétique lors du problème inverse. Le temps de calcul des 2003 problèmes inverses a été de 2 heures 10 minutes sur un pentiumIII à 450 MHz. Un autre test a également été réalisé en utilisant une simple interpolation trilinéaire pour le problème inverse, le temps de résolution des 2003 problèmes inverses est alors de 42 minutes. On note un gain de temps inférieur (facteur 3 seulement) à ce que pouvait espérer la figure J.2, le temps de résolution du problème inverse n'étant pas dans ce cas prépondérant par rapport au temps de chargement des grilles.

### 2.2 Résultats

Les erreurs de localisation obtenues sont représentées sur la figure J.6. On note que les erreurs sont sensiblement les mêmes pour l'interpolation trilinéaire et l'interpolation spline3D. L'écart avec le résultat obtenu par la méthode réaliste classique est de l'ordre de 1 mm pour les faibles profondeurs, mais tombe en dessous de 0.5 mm lorsque la profondeur dépasse 1 cm. Les cartes des erreurs de localisation en fonction de la profondeur sont représentées sur la figure J.7, la table J.2 fournit un récapitulatif de ces résultats.
<figure><img src="image_113.png" /><figcaption>FIG. J.6: Erreur de localisation du dipôle en fonction de la profondeur.</figcaption></figure>

Cette fois encore, l'utilisation de l'interpolation spline ne semble pas donner des résultats
meilleurs que l'interpolation trilinéaire.

<table><caption>TAB. J.2: Pourcentage des 2003 dipôles localisés avec moins de seuil mm selon le type de problème inverse utilisé</caption><thead><tr><td>seuil</td><td>grille,spline</td><td>grille,trilinéaire</td><td>BEM</td></tr></thead><tbody><tr><td>3 mm</td><td>94.7 %</td><td>93.8 %</td><td>97.0 %</td></tr><tr><td>6 mm</td><td>98.8 %</td><td>98.4 %</td><td>99.3 %</td></tr><tr><td>12 mm</td><td>99.6 %</td><td>99.6 %</td><td>99.7 %</td></tr></tbody></table>
<figure><img src="image_114.png" /><figcaption>FIG. J.7: Erreur de localisation du dipôle en fonction de sa position originale projetée sur l'IRM. La couleur est attribuée en fonction de l'erreur de localisation obtenue pour chaque dipôle : un point rouge représente un dipôle localisé avec une erreur inférieure à 3 mm ; jaune, 6 mm ; vert, 12 mm ; et bleu au delà.</figcaption></figure>
## 2.3 Influence du pas de la grille

Dans la section précédente, le problème inverse a été résolu avec un maillage plus grossier que celui utilisé pour créer les données simulées. Les erreurs obtenues pouvaient donc être dues à cette approximation géométrique plutôt qu'à l'utilisation des grilles. Pour avoir une meilleure idée du biais intrinsèque engendré par l'utilisation des grilles sur la localisation du dipôle, le maillage de référence a été utilisé dans l'étape de précalculs sur les grilles. L'influence du pas de la grille sur la minimisation a également été évalué.

### 2.3.1 Méthodes

Deux grilles ont été créées à partir du maillage de référence 3L1500 : l'une avec un pas de 10 mm et une troncature à 2 mm du maillage, et l'autre avec un pas de 5 mm et une troncature à 0.1 mm du maillage. Le calcul des champs engendrés pour ces grilles sur un pentiumIII-450 était de 16h35 pour la grille de pas 5 mm et de 1h25 pour la grille de pas 10 mm. Le temps de calcul des 2003 problèmes inverse dans le cas d'une interpolation spline était de 3h15. Les problèmes inverse sont également résolus avec la méthode par intégrale de surface classique pour le maillage 3L1500, afin de contrôler la part des erreurs dûes à la troncature de la grille. Le temps de résolution des problèmes inverses par l'approche classique a été de 46h50.

### 2.3.2 Résultats

Les résultats obtenus pour ces deux densités de grille avec les interpolations linéaire et spline sont donnés figure J.8 et table J.3. On constate que le pas de la grille a beaucoup plus d'influence sur l'erreur résiduelle que le type d'interpolation utilisé.

<table><caption>TAB. J.3: Pourcentage des 2003 dipôles localisés avec moins de seuil mm selon le type de problème inverse utilisé</caption><thead><tr><th>seuil</th><th>grille 10 mm, spline</th><th>grille 10 mm, trilinéaire</th><th>grille 5 mm, spline</th><th>grille 5 mm, trilinéaire</th></tr></thead><tbody><tr><td>1 mm</td><td>77.8 %</td><td>77.5 %</td><td>91.9 %</td><td>93.3 %</td></tr><tr><td>3 mm</td><td>94.7 %</td><td>95.0 %</td><td>98.7 %</td><td>99.0 %</td></tr><tr><td>6 mm</td><td>98.6 %</td><td>98.7 %</td><td>99.6 %</td><td>99.6 %</td></tr><tr><td>12 mm</td><td>99.7 %</td><td>99.7 %</td><td>99.7 %</td><td>99.8 %</td></tr></tbody></table>
<figure><img src="image_115.png" /><figcaption>FIG. J.8: Erreur résiduelle obtenue pour les différents pas de grille.</figcaption></figure>
<figure><img src="image_116.png" /><figcaption>FIG. J.9: Intervalles de confiance estimés pour la localisation d'un dipôle. Pour estimer les intervalles de confiance, le champ magnétique engendré par le dipôle équivalent trouvé est bruité, et pour chacune des 400 réalisations de bruit, le problème inverse est à nouveau résolu.</figcaption></figure>

## 2.4 Problème rencontré lors du calcul des intervalles de confiance

L'un des intérêts de l'approche proposée est une accélération des calculs de résolution du problème inverse permettant en pratique de calculer des intervalles de confiance en modifiant les données du problème (bruit ajouté aux données) pour un grand nombre de réalisations. On notera toutefois que l'approche proposée ne permet pas d'évaluer un intervalle de confiance en bruitant les positions des voies, puisqu'il faudrait à chaque réalisation de bruit refaire les calculs préliminaires.

Nous avons donc cherché à évaluer l'approche proposée dans le cadre du calcul des intervalles de confiance réalisés en bruitant uniquement les données de potentiel et de champ magnétique. Les intervalles de confiance auxquels nous ferons référence dans la suite de cette annexe ont été calculés en ajoutant un bruit gaussien dont le rapport signal sur bruit est de 5. Pour chaque intervalle de confiance, 400 réalisations de bruit ont été effectuées. L'intervalle de confiance est ensuite représenté soit par la totalité des points correspondant aux localisations de sources obtenues pour chaque réalisation de bruit, soit par l'ellipsoïde approchant au mieux l'espace occupé par ces points.