# Méthodes de localisation des générateurs de l'activité électrique cérébrale à partir de signaux électro- et magnéto-encéphalographiques.

Anne Eugénie, Jeanne Cheylus Née Crouzeix

## ► To cite this version:

Anne Eugénie, Jeanne Cheylus Née Crouzeix. Méthodes de localisation des générateurs de l'activité électrique cérébrale à partir de signaux électro- et magnéto-encéphalographiques. Traitement du signal et de l'image [eess.SP]. Institut national des sciences appliquées de Lyon, 2001. Français. NNT: tel-03826298

HAL Id: tel-03826298

https://theses.hal.science/tel-03826298

Submitted on 24 Oct 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
THÈSE

présentée

DEVANT L'INSTITUT NATIONAL DES SCIENCES APPLIQUÉES DE LYON

pour obtenir

LE GRADE DE DOCTEUR

FORMATION DOCTORALE : IMAGES ET SYSTÈMES

ÉCOLE DOCTORALE : ÉLECTRONIQUE, ÉLECTROTECHNIQUE, AUTOMATIQUE

par

Anne CROUZEIX épouse CHEYLUS

Ingénieur de l'École Centrale de Lyon

Méthodes de localisation des générateurs de l'activité électrique cérébrale à partir de signaux électro- et magnéto-encéphalographiques.

Soutenue le 8 octobre 2001 devant la Commission d'examen

Jury MM.

Olivier BERTRAND

Directeur de recherche

INSERM U280

directeur de thèse

Line GARNERO

Directeur de recherche

LENA CNRS

rapporteur

Isabelle MAGNIN

Directeur de recherche

INSA LYON

examenateur

Michael SCHERG

Professeur

Université de HEIDELBERG

rapporteur

Blaise YVERT

Chargé de recherche

INSERM U280

examenateur

Cette thèse a été préparée au laboratoire "Processus Mentaux et Activation Cérébrale", INSERM U280, LYON.
<figure><img src="image_2.png" /></figure>
Remerciements

Tout d'abord, je voudrais remercier chaleureusement Olivier BERTRAND, sans qui je
n'aurais jamais eu l'idée de faire une thèse sur ce sujet. Il m'a initiée non seulement à la
modélisation de l'activité cérébrale en MEG et en EEG, mais également à de nombreux
aspects du travail de chercheur. Par son suivi en tant que directeur de thèse et ses
encouragements, il a su conserver ce projet dans la bonne voie.

Cet encadrement a été largement accompagné par Jacques PERNIER et Blaise YVERT,
que je voudrais également remercier ici. Blaise a suivi avec attention toute l'évolution de ce
travail, et ses commentaires, ses relectures et son enthousiasme m'ont été très précieux. Le
suivi et les relectures de Jacques PERNIER m'ont également permis de renforcer la rigueur de
ce travail de thèse, ses conseils ont été particulièrement utiles pour les aspects théoriques de
cette thèse.

Je voudrais également remercier Line GARNERO, pour m'avoir fait parvenir les
enregistrements simultanés de MEG et d'EEG qui m'ont permis de générer un bruit réaliste et
corrélé spatialement pour valider mon algorithme de couplage au chapitre 4.

Je souhaite également remercier Jean-François ECHALLIER et Pierre-Emmanuel
AGUERA pour leur précieuse aide technique et informatique. Merci pour leur patience, leur
calme et leur sens de l'humour lorsque réseau, imprimantes et ordinateurs adoptaient un
comportement chaotique. Merci également pour leur ouverture vis à vis de mon goût pour les
logiciels libres.

Je tiens aussi à remercier Martine THEILLERE pour avoir su accomplir toutes sortes de
démarches administratives avec une bonne humeur communicative.

Je remercie l'ensemble du laboratoire de l'unité 280 pour son accueil chaleureux et tous
les bons moments que j'y ai passés, que ce soit dans le cadre de grandes discussions ayant trait
à la recherche ou à refaire le monde ou encore dans le cadre de sympathiques vendanges.

Je remercie également mon mari pour son soutien constant.
<figure><img src="image_3.png" /></figure>
# Avant propos

Ce travail de thèse porte sur les moyens d'améliorer les méthodes de localisation des sources neuronales de l'activité électrique cérébrale telle qu'elle peut être mesurée à l'extérieur de la tête par électroencéphalographie (EEG) ou par magnétoencéphalographie (MEG). L'évaluation de l'apport de l'utilisation de modèles réalistes avec des données magnétoencéphalographiques d'une part, et l'étude du couplage des données électro- et magnéto-encéphalographiques d'autre part sont les deux objectifs de cette thèse.

Dans un premier temps, au chapitre 1, le positionnement et les enjeux de cette recherche sont discutés, les équations fondamentales du problème sont rappelées, puis une revue bibliographique permet de préciser l'état de l'art sur l'évaluation de différentes approches de localisation des sources de l'activité électrique cérébrale.

Le chapitre 2 présente une première tentative d'amélioration de la localisation des sources à partir de données électriques en utilisant une fonction d'interpolation plus précise pour les calculs numériques par la méthodes des intégrales de surface.

Au chapitre 3, une évaluation de la méthode de localisation des sources à partir de données magnétiques suivant que l'on utilise un modèle sphérique ou plus réaliste de la tête est réalisée.

L'apport de l'utilisation simultanée des deux modalités électrique et magnétique est recherché au chapitre 4. Plusieurs méthodes de couplage sont évaluées et comparées à l'utilisation séparée des données de chaque modalité.

Enfin, les conclusions et perspectives de ce travail de thèse font l'objet du chapitre 5.
<figure><img src="image_4.png" /></figure>
# Table des matières

Avant propos i

1 Introduction 3

1.1 Généralités 4

1.1.1 La localisation des sources de l'activité électrique cérébrale 4

1.1.2 Activité électrique au niveau cellulaire. 7

1.1.3 Modélisation de l'activité électrique cérébrale . . . . . . . . . . . . . 10

1.2 Formalisation mathématique du problème . . . . . . . . . . . . . . . . . . . . 13

1.2.1 Équations générales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.2.2 Caractéristiques des champs de potentiel et des champs magnétiques . . 15

1.2.3 La méthode des intégrales de surface . . . . . . . . . . . . . . . . . . 16

1.2.4 Méthodes de résolution du problème inverse . . . . . . . . . . . . . . . . 20

1.2.5 Mise en œuvre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

1.3 Synthèse bibliographique de l'évaluation des modèles . . . . . . . . . . . . . 30

1.3.1 Apport des modèles réalistes par rapport aux modèles sphériques . . . . 30

1.3.2 Apport de l'interpolation linéaire du potentiel dans la méthode des intégrales du surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

1.3.3 Apport de l'affinage local des maillages . . . . . . . . . . . . . . . . . . . 36

1.3.4 Apport du couplage MEG/EEG . . . . . . . . . . . . . . . . . . . . . . . 37

1.4 Plan de travail . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

2 Amélioration de la méthode BEM en EEG par une interpolation plus précise 41

2.1 L'interpolation linéaire du potentiel . . . . . . . . . . . . . . . . . . . . . . . . 42

2.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

2.1.2 Méthode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

2.1.3 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

2.1.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

2.2 Amélioration de l'interpolation géométrique . . . . . . . . . . . . . . . . . . . . 54

2.3 Conclusions et perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3 Apport des modèles réalistes en MEG 57

3.1 Introduction 58

3.1.1 Revue bibliographique 58

3.1.2 plan de travail 62

3.2 Évaluation de l'apport des modèles réalistes pour la résolution du problème inverse 63

3.2.1 Méthode 63

3.2.2 Résultats 70

3.2.3 Discussion 74

3.3 Effet de l'orientation du dipôle 80

3.3.1 Méthode 81

3.3.2 Topographie du champ magnétique engendré par un dipôle radial dans un modèle réaliste 82

3.3.3 Effet de l'orientation des sources sur l'amplitude du signal mesuré 86

3.3.4 Effet de l'orientation des sources sur les erreurs de localisation et d'orientation. 88

3.4 Conclusions et perspectives 88

4 Évaluation de plusieurs méthodes de couplage MEG/EEG 91

4.1 Introduction 92

4.2 Objectifs 95

4.3 Paramètre d'ajustement de la conductivité 96

4.4 Modèle de bruit 99

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

4.5 Modèle utilisé pour la génération des signaux bruités 100

4.6 Évaluation de l'apport du couplage 100

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

4.6.2 Montages utilisés 101

4.6.3 Critères à minimiser 101

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

4.7 Première configuration : sources à plusieurs profondeurs 106

4.7.1 Résultats 106

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

4.8.1 Résultats 112

4.9 Troisième configuration : une source dans chaque hémisphère 117

4.9.1 Résultats 119

4.10 Quatrième configuration : deux sources dans le même hémisphère 119
<table><tr><td>4.10.1 Résultats</td><td>119</td></tr><tr><td>4.11 Conclusions</td><td>122</td></tr></table>

<table><tr><td>5 Conclusions et perspectives</td><td>123</td></tr><tr><td>5.1 Conclusions</td><td>124</td></tr><tr><td>5.2 Perspectives</td><td>125</td></tr></table>

<table><tr><td>Annexes</td><td>129</td></tr></table>

<table><tr><td>A Le potentiel et le champ magnétique dans le modèle sphérique</td><td>129</td></tr><tr><td>1 Le potentiel électrique [Ary et al., 1981]</td><td>129</td></tr><tr><td>2 Le champ magnétique [Sarvas, 1987]</td><td>130</td></tr></table>

<table><tr><td>B Équations intégrales du potentiel et du champ magnétique</td><td>133</td></tr><tr><td>1 Équation intégrale du potentiel</td><td>133</td></tr><tr><td>1.1 Résolution de Δu = f avec u ∈ D</td><td>133</td></tr><tr><td>1.2 Application à l'équation de Poisson [Geselowitz, 1967]</td><td>133</td></tr><tr><td>2 Équation intégrale du champ magnétique [Geselowitz, 1970]</td><td>135</td></tr></table>

<table><tr><td>C Calcul des matrices dans la méthode BEM</td><td>137</td></tr><tr><td>1 Calcul des matrices BEM pour le potentiel électrique</td><td>137</td></tr><tr><td>1.1 L'approche centre de masse [Hämäläinen, 1987]</td><td>137</td></tr><tr><td>1.2 L'approche interpolation linéaire [De Munck, 1992]</td><td>138</td></tr><tr><td>2 Calcul des matrices BEM pour le champ magnétique</td><td>142</td></tr><tr><td>2.1 Approche centre de masse</td><td>143</td></tr><tr><td>2.2 Approche interpolation linéaire [Ferguson et al., 1994]</td><td>143</td></tr></table>

<table><tr><td>D Influence de la conductivité en EEG et en MEG</td><td>147</td></tr><tr><td>1 Effet d'un facteur multiplicatif affectant toutes les conductivités en EEG</td><td>147</td></tr><tr><td>2 Effet d'un facteur multiplicatif affectant toutes les conductivités en EEG</td><td>148</td></tr></table>

<table><tr><td>E Forme matricielle de l'approche problème isolé</td><td>149</td></tr><tr><td>1 Introduction</td><td>149</td></tr><tr><td>2 Problème à p couches.</td><td>150</td></tr><tr><td>3 Passage au problème à n couches</td><td>150</td></tr></table>

<table><tr><td>F Calcul des matrices dans l'approche LINT BEM avec des triangles courbes</td><td>153</td></tr><tr><td>1 Approche LINT sur un triangle plat</td><td>153</td></tr><tr><td>2 Approche LINT sur un triangle courbe</td><td>153</td></tr></table>
G Article paru dans : Clinical Neurophysiology 110 (1999) 2176-2188 157

<table><tr><td>1</td><td>Introduction</td><td>158</td></tr><tr><td>2</td><td>Methods</td><td>162</td></tr><tr><td>2.1</td><td>Mesh generation</td><td>162</td></tr><tr><td>2.2</td><td>Forward model</td><td>162</td></tr><tr><td>2.3</td><td>MEG sensor configuration</td><td>164</td></tr><tr><td>2.4</td><td>Choice of the dipole sites</td><td>164</td></tr><tr><td>2.5</td><td>Inverse models</td><td>166</td></tr><tr><td>2.6</td><td>The inverse procedure</td><td>167</td></tr><tr><td>2.7</td><td>Error criteria</td><td>167</td></tr><tr><td>3</td><td>Results</td><td>172</td></tr><tr><td>3.1</td><td>Bias due to the spherical approximation</td><td>172</td></tr><tr><td>3.2</td><td>Errors obtained with realistic models</td><td>175</td></tr><tr><td>4</td><td>Discussion</td><td>176</td></tr><tr><td>4.1</td><td>Localization bias due to spherical model</td><td>177</td></tr><tr><td>4.2</td><td>Influence of the number of layers with realistic models.</td><td>178</td></tr><tr><td>4.3</td><td>Orientation and moment amplitude errors</td><td>178</td></tr><tr><td>5</td><td>References</td><td>179</td></tr></table>

H Méthode de Choleski 183

<table><tr><td>1</td><td>Cas général</td><td>183</td></tr><tr><td>2</td><td>Cas des matrices symétriques positives non définies</td><td>185</td></tr></table>

I Génération d'un bruit réaliste 191

<table><tr><td>1</td><td>Traitement d'un enregistrement couplé spontané pour générer un bruit réaliste</td><td>191</td></tr><tr><td>1.1</td><td>Recalage des électrodes et des canaux MEG sur un modèle réaliste</td><td>191</td></tr><tr><td>1.2</td><td>Traitement des signaux</td><td>193</td></tr></table>

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

<table><tr><td>1</td><td>Présentation et évaluation de la méthode dans le cas d'une géométrie sphérique</td><td>199</td></tr><tr><td>1.1</td><td>Introduction</td><td>199</td></tr><tr><td>1.2</td><td>Méthodes</td><td>200</td></tr><tr><td>1.3</td><td>Résultats</td><td>203</td></tr><tr><td>2</td><td>Validation préliminaire dans le cas d'une géométrie réaliste en MEG</td><td>209</td></tr><tr><td>2.1</td><td>Méthodes</td><td>209</td></tr><tr><td>2.2</td><td>Résultats</td><td>209</td></tr><tr><td>2.3</td><td>Influence du pas de la grille</td><td>212</td></tr><tr><td>2.4</td><td>Problème rencontré lors du calcul des intervalles de confiance</td><td>214</td></tr></table>
<table><tr><td>2.5</td><td>Recherche de la cause des accumulations linéaires des intervalles de confiance</td><td>215</td></tr><tr><td>2.6</td><td>Validité des grilles pour l'utilisation des intervalles de confiance</td><td>216</td></tr><tr><td>2.7</td><td>Conclusion</td><td>220</td></tr><tr><td>K</td><td>Liste des publications issues de cette thèse</td><td>221</td></tr><tr><td>1</td><td>Publications</td><td>221</td></tr><tr><td>2</td><td>Résumés publiés</td><td>221</td></tr><tr><td>3</td><td>Communications à des congrès</td><td>222</td></tr><tr><td>4</td><td>Rapports techniques</td><td>223</td></tr><tr><td></td><td>Bibliographie</td><td>225</td></tr><tr><td></td><td>Index</td><td>236</td></tr></table>
<figure><img src="image_5.png" /></figure>
# Liste des tableaux

2.1 Maillages utilisés pour le problème direct. GD = densité globale (triangles/cm²). LD = densité locale (triangles/cm²). N = nombre de triangles. . 45

2.2 Maillages utilisés pour le problème inverse. GD = densité globale (triangles/cm²). LD = densité locale (triangles/cm²). N = nombre de triangles. . 45

3.1 Erreurs de localisation en modèle sphérique selon la littérature . . . . . . . . . 61

3.2 Maillages utilisés pour le problème direct et le problème inverse . . . . . . . . 65

G.1 Localization errors found in the literature . . . . . . . . . . . . . . . . . . . . . 161

G.2 Characteristics of the models used for the forward and inverse problem. . . . . 162

J.1 Temps de calcul des champs engendrés . . . . . . . . . . . . . . . . . . . . . . . . 204

J.2 Pourcentage des 2003 dipôles localisés avec moins de seuil mm . . . . . . . . 210

J.3 Pourcentage des 2003 dipôles localisés avec moins de seuil mm . . . . . . . . 212
<figure><img src="image_6.png" /></figure>
Table des figures

<table>
  <tbody>
    <tr>
      <td>1.1</td>
      <td>Description d'un neurone pyramidal et d'un neurone étoilé.</td>
      <td>9</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td>Schéma de propagation du potentiel d'action</td>
      <td>9</td>
    </tr>
    <tr>
      <td>1.3</td>
      <td>Exemple d'utilisation du modèle sphérique</td>
      <td>12</td>
    </tr>
    <tr>
      <td>1.4</td>
      <td>Exemple de maillage réaliste de la tête</td>
      <td>12</td>
    </tr>
    <tr>
      <td>1.5</td>
      <td>Algorithme de Marquardt [1963] modifié par Nash [1979]</td>
      <td>25</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Exemple de maillage de sphère localement affiné sur un angle d'ouverture de 40°</td>
      <td>44</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Position des électrodes sur le scalp.</td>
      <td>46</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Positions des axes temporal et pariétal</td>
      <td>47</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Comparaisons maillages homothétiques ou non pour l'approche COM</td>
      <td>48</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>Comparaison entre LINT et COM pour le problème direct</td>
      <td>49</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>MAG et RDM en fonction du modèle</td>
      <td>50</td>
    </tr>
    <tr>
      <td>2.7</td>
      <td>MAG et RDM en fonction de la taille de la matrice</td>
      <td>50</td>
    </tr>
    <tr>
      <td>2.8</td>
      <td>Érreur de localisation en fonction du modèle</td>
      <td>51</td>
    </tr>
    <tr>
      <td>2.9</td>
      <td>Érreur de localisation en fonction du nombre d'électrodes</td>
      <td>52</td>
    </tr>
    <tr>
      <td>2.10</td>
      <td>Érreur de localisation en fonction de la taille de la matrice</td>
      <td>52</td>
    </tr>
    <tr>
      <td>2.11</td>
      <td>maillages à triangles plats et à triangles courbes</td>
      <td>54</td>
    </tr>
    <tr>
      <td>2.12</td>
      <td>Érreur de localisation pour les triangles plats et courbes dans l'approche LINT.</td>
      <td>55</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Procédure de simulation</td>
      <td>64</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Positions des 2003 dipôles</td>
      <td>66</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>RMS des 2003 dipôles</td>
      <td>66</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Tangentialités des 2003 dipôles</td>
      <td>67</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Procédure adoptée pour choisir une sphère locale</td>
      <td>68</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Érreurs pour les modèles sphériques en MEG</td>
      <td>71</td>
    </tr>
    <tr>
      <td>3.7</td>
      <td>Érreur en fonction de la tangentialité</td>
      <td>73</td>
    </tr>
    <tr>
      <td>3.8</td>
      <td>Érreur de localisation (vues sur IRM)</td>
      <td>75</td>
    </tr>
    <tr>
      <td>3.9</td>
      <td>Érreur d'orientation (vues sur IRM)</td>
      <td>76</td>
    </tr>
    <tr>
      <td>3.10</td>
      <td>Érreurs selon la position verticale en MEG et en EEG</td>
      <td>77</td>
    </tr>
    <tr>
      <td>3.11</td>
      <td>Dipôles de la simulation : une position, 180 orientations.</td>
      <td>81</td>
    </tr>
    <tr>
      <td>3.12</td>
      <td>Topographies engendrés par des dipôles radiaux en MEG</td>
      <td>83</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>3.13</td>
      <td>Évolution de la topographie en fonction de l'angle du dipôle : de 30° en 30°</td>
      <td>84</td>
    </tr>
    <tr>
      <td>3.14</td>
      <td>Évolution de la topographie en fonction de l'angle du dipôle : de 5° en 5°</td>
      <td>85</td>
    </tr>
    <tr>
      <td>3.15</td>
      <td>Évolution du RMS en fonction de l'angle du dipôle</td>
      <td>87</td>
    </tr>
    <tr>
      <td>3.16</td>
      <td>Exemples de rapports de RMS radial/tangentiel sur un axe temporal.</td>
      <td>87</td>
    </tr>
    <tr>
      <td>3.17</td>
      <td>Érreur de localisation en fonction de l'orientation du dipôle</td>
      <td>89</td>
    </tr>
    <tr>
      <td>3.18</td>
      <td>Érreur d'orientation en fonction de l'orientation du dipôle</td>
      <td>89</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Paramètre d'ajustement de la conductivité</td>
      <td>98</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Représentation graphique de l'intercorrélation normalisée entre les voies.</td>
      <td>99</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Modèle utilisé pour la méthode des intégrales de surface.</td>
      <td>100</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Positions des électrodes et des voies MEG.</td>
      <td>102</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Première configuration des sources</td>
      <td>106</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>Rapports signal sur bruit obtenus pour la première configuration de sources.</td>
      <td>107</td>
    </tr>
    <tr>
      <td>4.7</td>
      <td>4 réalisations de signal bruité</td>
      <td>107</td>
    </tr>
    <tr>
      <td>4.8</td>
      <td>Érreurs pour la première configuration de source et bruit1</td>
      <td>108</td>
    </tr>
    <tr>
      <td>4.9</td>
      <td>Érreurs pour la première configuration de sources et bruit2</td>
      <td>109</td>
    </tr>
    <tr>
      <td>4.10</td>
      <td>Deuxième configuration des sources</td>
      <td>111</td>
    </tr>
    <tr>
      <td>4.11</td>
      <td>Rapports signal sur bruit obtenus pour la deuxième configuration de sources.</td>
      <td>112</td>
    </tr>
    <tr>
      <td>4.12</td>
      <td>Érreurs pour la deuxième configuration de sources</td>
      <td>113</td>
    </tr>
    <tr>
      <td>4.13</td>
      <td>Érreurs de localisation et d'orientation en MEG réaliste ou sphérique bruitée</td>
      <td>116</td>
    </tr>
    <tr>
      <td>4.14</td>
      <td>Troisième configuration des sources</td>
      <td>117</td>
    </tr>
    <tr>
      <td>4.15</td>
      <td>Rapports signal sur bruit obtenus pour la troisième configuration de sources.</td>
      <td>117</td>
    </tr>
    <tr>
      <td>4.16</td>
      <td>Érreurs pour la troisième configuration de sources</td>
      <td>118</td>
    </tr>
    <tr>
      <td>4.17</td>
      <td>Quatrième configuration des sources</td>
      <td>120</td>
    </tr>
    <tr>
      <td>4.18</td>
      <td>Rapports signal sur bruit obtenus pour la quatrième configuration de sources.</td>
      <td>120</td>
    </tr>
    <tr>
      <td>4.19</td>
      <td>Érreurs pour la quatrième configuration de sources</td>
      <td>121</td>
    </tr>
    <tr>
      <td>A.1</td>
      <td>Description du modèle sphérique en EEG</td>
      <td>129</td>
    </tr>
    <tr>
      <td>G.1</td>
      <td>Simulation procedure</td>
      <td>163</td>
    </tr>
    <tr>
      <td>G.2</td>
      <td>Dipole positions, RMS and tangentiality</td>
      <td>165</td>
    </tr>
    <tr>
      <td>G.3</td>
      <td>Errors vs slice number and tangentiality, spherical models</td>
      <td>169</td>
    </tr>
    <tr>
      <td>G.4</td>
      <td>Errors vs tangentiality, realistic models</td>
      <td>170</td>
    </tr>
    <tr>
      <td>G.5</td>
      <td>Localization error for BRAINLOCSPH and realistic models.</td>
      <td>171</td>
    </tr>
    <tr>
      <td>G.6</td>
      <td>Orientation error for spherical and realistic models.</td>
      <td>173</td>
    </tr>
    <tr>
      <td>G.7</td>
      <td>Errors vs slice number for MEG and EEG</td>
      <td>174</td>
    </tr>
    <tr>
      <td>I.1</td>
      <td>Positions des voies EEG.</td>
      <td>192</td>
    </tr>
    <tr>
      <td>I.2</td>
      <td>Détermination du repère 10-20 à partir des positions des électrodes.</td>
      <td>192</td>
    </tr>
    <tr>
      <td>I.3</td>
      <td>Passage du repère 10-20 au repère tête de l'unité 280.</td>
      <td>193</td>
    </tr>
    <tr>
      <td>I.4</td>
      <td>Positions des voies MEG (en rouge) et EEG (en bleu) recalées sur notre modèle.</td>
      <td>194</td>
    </tr>
  </tbody>
</table>
<table><tr><td>I.5</td><td>Spectres d'une voie MEG (à gauche) et d'une voie d'EEG (à droite).</td><td>194</td></tr><tr><td>I.6</td><td>Cartes des écarts types par voie des réalisations de bruit</td><td>195</td></tr><tr><td>I.7</td><td>Cartes de la première composante principale</td><td>196</td></tr><tr><td>I.8</td><td>Topographies de quatre réalisations de bruit.</td><td>197</td></tr><tr><td>I.9</td><td>Représentation graphique de l'intercorrélation normalisée entre les voies.</td><td>197</td></tr><tr><td>I.10</td><td>Positions des électrodes dans les configurations à 29, 64, 151 ou 256 électrodes.</td><td>198</td></tr><tr><td>J.1</td><td>Positions des voies, grille, et positions des sources</td><td>202</td></tr><tr><td>J.2</td><td>Temps de calcul pour les différentes approches.</td><td>204</td></tr><tr><td>J.3</td><td>Erreurs intrinsèques à la méthode</td><td>206</td></tr><tr><td>J.4</td><td>Grilles non tronquées : erreurs obtenues</td><td>207</td></tr><tr><td>J.5</td><td>Grilles tronquées : erreurs obtenues</td><td>208</td></tr><tr><td>J.6</td><td>Errore de localisation du dipôle en fonction de la profondeur.</td><td>210</td></tr><tr><td>J.7</td><td>Erreurs de localisation des dipôles</td><td>211</td></tr><tr><td>J.8</td><td>Errour résiduelle obtenue pour les différents pas de grille.</td><td>213</td></tr><tr><td>J.9</td><td>Intervalles de confiance estimés pour la localisation d'un dipôle.</td><td>214</td></tr><tr><td>J.10</td><td>Errour résiduelle comparée pour les différentes interpolations.</td><td>217</td></tr><tr><td>J.11</td><td>Intervalles de confiance pour un dipôle de coordonnées (9,9,69,1,-1,0)</td><td>218</td></tr><tr><td>J.12</td><td>Intervalles de confiance pour un dipôle de coordonnées (3,-11,98,0,0,1)</td><td>219</td></tr></table>
<figure><img src="image_7.png" /></figure>
Il réduisit l'écartement de leurs deux corps par le moyen d'un raccourcissement du biceps droit transmis du cerveau, le long d'une paire de nerfs crâniens choisie judicieusement.

Boris Vian, *l'écume des jours*.
<figure><img src="image_8.png" /></figure>
# Chapitre 1

## Introduction
# 1.1 Généralités

## 1.1.1 La localisation des sources de l'activité électrique cérébrale

### 1.1.1.1 Électroencéphalogramme

L'électroencéphalogramme (EEG) est l'enregistrement de différences de potentiels électriques effectués sur la surface de la tête.

D'après Changeux [1983], l'activité électrique cérébrale fut enregistrée pour la première fois par un jeune assistant de physiologie de la Royal Infirmary de Liverpool nommé Caton en 1875. En ouvrant la boîte crânienne d'un singe, il observe la présence de courants électriques attestés par les oscillations d'un galvanomètre. Lorsque la rétine est soumise à un stimulus lumineux, il observe une variation plus importante du courant électrique localisée dans la région occipitale du cortex. En 1929, le neuropsychiatre Allemand Hans Berger met au point une technique d'enregistrement permettant de mesurer l'activité électrique cérébrale sans ouvrir la boîte crânienne. Il utilise des électrodes posées sur le scalp et reliées à des amplificateurs pour mesurer des oscillations électriques de quelques dizaines de microvolts.

Cette nouvelle technique, l'électroencéphalographie, donnera lieu à de nombreuses études de l'activité électrique cérébrale. On observera tout d'abord des oscillations de forte amplitude et ayant une fréquence de l'ordre de 10 Hz, et dont l'amplitude décroît lorsque l'on focalise son attention. Ce sont les oscillations alpha. L'analyse fréquentielle de l'EEG permet de mettre en évidence d'autres rythmes (delta : inférieur à 4 Hz, thêta : entre 4 et 8 Hz, bêta : entre 8 et 20 Hz). L'apparition de ces rythmes, ainsi que d'autres signaux, interviennent dans la classification des différents stades d'éveil et de sommeil. Une autre façon d'analyser l'EEG consiste à stimuler de manière répétée un sujet et à faire la moyenne des signaux électriques mesurés sur chaque électrode dans une certaine tranche de temps entourant le stimulus. Ce processus permet d'éliminer en grande partie les signaux électriques liés à l'activité spontanée du cerveau, et d'extraire uniquement la réponse liée au stimulus. C'est ce que l'on appelle les potentiels évoqués.

### 1.1.1.2 Magnétoencéphalogramme

Le magnétoencéphalogramme (MEG) est la représentation graphique des variations de champ magnétique mesurées au voisinage de la tête. Du fait de la très faible intensité des champs magnétiques résultant de l'activité électrique cérébrale, les premières expériences de MEG furent beaucoup plus tardives que l'EEG. Cohen [1968] parvint tout de même à mesurer avec une spire classique l'activité ayant la plus forte amplitude (les oscillations alpha). Ce n'est cependant qu'avec l'avènement des supraconducteurs que l'on parvint à construire un
capteur, constitué de spires reliées à un SQUID (Superconducting Quantum Interference Device), et pouvant mesurer les champs magnétiques cérébraux de l'ordre du femto-Tesla ($1 \text{ fT} = 10^{-15} \text{ T}$). A l'aide de cette technologie, Cohen [1972] parvint à faire le premier enregistrement de MEG.

### 1.1.1.3 Localisation des sources de l'activité électrique cérébrale

Toujours selon Changeux [1983], le lien entre la perte de certaines facultés cognitives telles que le langage, et l'atteinte de certaines zones du cerveau avait été fait dès l'antiquité égyptienne. Toutefois, le débat persistera pendant une très longue période pour savoir si le siège de l'âme se trouve dans le cœur ou dans le cerveau, et ce n'est qu'au $XIX^{\text{e}}$ siècle que la théorie de la localisation cérébrale des fonctions mentales émergera. Vers 1800, l'Allemand Franz Joseph Gall découpe le cerveau en plusieurs zones, auxquelles il associe des facultés mentales. Si les localisations et le découpage des fonctions mentales étaient plutôt fantaisistes - ainsi pensait-il que la palpation des bosses du crâne permettait d'évaluer la vertu, la générosité ou l'avarice d'une personne - l'idée de la localisation de zones cérébrales associées à une faculté mentale a créé un véritable courant dans la recherche neurologique. En 1861, Paul Broca établit d'une manière beaucoup plus scientifique un lien entre une région du cerveau et une fonction cérébrale. Cette découverte lui est venue lors de l'autopsie du cerveau d'un de ses patients, surnommé "Tan", car il souffrait d'un trouble du langage articulé et ne pouvait prononcer que ce phonème. Le cerveau de ce patient montrait une lésion dans une région arrière du lobe frontal aujourd'hui appelée aire de Broca. Par la suite, de nouvelles techniques d'imagerie cérébrale confirmeront ce lien entre l'aire de Broca et les mécanismes cérébraux impliqués dans la production du langage articulé. Nous avons vu que Caton, dans ses expériences, a également mis en évidence une localisation des aires cérébrales impliquées dans la perception d'un stimulus lumineux. L'électrophysiologie à crâne ouvert, la stimulation électrique chez l'animal, l'autopsie et l'étude des lésions cérébrales chez l'homme ont pendant de longues années été les seuls outils de la localisation d'aires fonctionnelles cérébrales.

Aujourd'hui, de nouvelles techniques sont apparues, permettant une localisation non traumatique in vivo des fonctions cérébrales chez l'homme. Ce sont la tomographie à émission de positons (TEP), l'imagerie par résonance magnétique fonctionnelle (IRMf), la tomographie par émission de simples photons (SPECT), ainsi que la MEG et l'EEG. La TEP consiste à injecter dans le sang un marqueur radioactif, puis à détecter tout autour de la tête les paires de photons gamma qui sont émis lorsqu'un positon s'annihile avec un électron. Un algorithme de reconstruction permet d'aboutir à une carte en trois dimensions des zones d'émission. La fonction physiologique étudiée dépend du type de molécule marquée. Par exemple, l'eau marquée par l'oxygène 15 permet d'estimer les variations de débit sanguin
local induites par une activité du cerveau. La résolution spatiale de cette technique est de l'ordre de 4 mm pour un enregistrement qui dure de 1 à 2 minutes. La SPECT fait également appel à des traceurs radioactifs, mais les réactions de désintégration ne donnent lieu qu'à l'émission d'un seul photon gamma. Il en résulte une localisation plus difficile des lieux de l'émission, et une résolution spatiale plus faible. L'IRM fonctionnelle utilise la fréquence de résonance du moment magnétique, ou spin, de noyaux atomiques comme l'hydrogène et un champ magnétique variable dans l'espace pour reconstruire des images en trois dimensions. Selon la technique utilisée, on obtient des images anatomiques très précises du cerveau, avec une résolution spatiale inférieure à 1 mm, ou des images fonctionnelles avec un temps d'acquisition beaucoup plus court, jusqu'au 50ᵉ de seconde, mais avec une résolution spatiale plus faible (3 mm), et correspondant généralement aux zones où le taux d'oxyhémoglobine augmente.

Par contraste, les techniques classiques de reconstruction des sources neuronales à partir de l'EEG, qui utilisent un modèle de la tête dont nous reparlerons à la section suivante (section 1.1.3.2), ne permettent d'obtenir qu'une résolution spatiale de l'ordre du centimètre, et pour la MEG, de l'ordre du demi-centimètre, dans l'hypothèse pas toujours vérifiée d'une activité cérébrale se modélisant bien par un simple dipôle de courant. La résolution temporelle de la MEG et de l'EEG, de l'ordre de la milliseconde, est bien supérieure à celle des autres techniques d'imagerie cérébrale. L'amélioration des techniques de reconstruction des sources de l'activité électrique cérébrale devrait donc permettre d'accéder à une information non seulement sur les lieux des activités cérébrales, mais plus globalement sur leur organisation spatio-temporelle.

#### 1.1.1.4 Perspectives d'application d'une meilleure localisation des sources

Une meilleure technique de localisation des sources de l'activité électrique cérébrale permettrait, nous l'avons vu, de fournir un nouvel outil pour la caractérisation des fonctions des structures cérébrales dans leur organisation tant spatiale que temporelle. Un tel outil présente un intérêt majeur pour la recherche sur les processus mentaux et cognitifs. En effet, les autres techniques d'investigation in vivo présentent une moins bonne résolution temporelle. L'EEG et la MEG présentent l'avantage d'être des techniques d'investigation non invasives et peu contraignantes pour le sujet car aucune injection n'est nécessaire, et les appareils ne sont pas bruyants. Pour ces raisons, de nombreuses études utilisent déjà l'une ou l'autre de ces modalités en proposant une localisation approximative des sources. Un exemple est l'étude de l'organisation tonotopique du cortex auditif [Romani et al., 1982; Pantev et al., 1988; Bertrand et al., 1988]. L'organisation somatotopique du cortex somato-sensoriel a également fait l'objet de ce type d'investigation, par exemple dans Buchner et al. [1994]. Mais
de nombreux autres sujets sont abordés, en particulier depuis qu’existent des logiciels commerciaux de localisation des sources à partir de l’EEG ou de la MEG (BESA®*, CURRY®*, ainsi que les logiciels fournis avec les installations de MEG).

De plus, il existe un certain nombre de cas cliniques où l’amélioration de la localisation des sources permettrait d’aider le praticien à faire un meilleur diagnostic de l’état de ses patients. Un exemple est le cas de patients souffrant d’une épilepsie rebelle à tout traitement médicamenteux. Il peut être envisagé dans ces cas là une ablation de la zone du cerveau dont le dérèglement est à l’origine de la crise épileptique, la zone épileptogène. Cette zone semble se caractériser par une activité électrique bien particulière, qui provoque l’activation des zones voisines du cerveau. Actuellement, la localisation de la zone épileptogène se fait à l’aide d’électrodes implantées dans le cerveau. Comme il s’agit souvent de zones cérébrales endommagées, la localisation d’une tumeur ou d’un traumatisme par IRM ou par scanner X fournit une information complémentaire. La précision de la localisation est dans ce cas cruciale, car plus l’ablation concerne une grande zone du cerveau, plus les risques de séquelles fonctionnelles sont importants. Une meilleure localisation du foyer épileptogène par MEG ou par EEG permettrait au minimum de mieux cibler la zone où implanter des électrodes intra-cérébrales pour une caractérisation plus précise de la zone incriminée, et dans le meilleur des cas, d’éviter d’avoir à implanter des électrodes intra-cérébrales. Cette application a été par exemple étudiée dans Roth et al. [1997]; Waberski et al. [1998]; Silva et al. [1999].

Un autre exemple d’application clinique pourrait être le suivi de la réorganisation fonctionnelle cérébrale dans le cas de patients ayant eu une lésion d’une zone du cerveau associée à une certaine activité mentale. On a observé, dans des cas similaires, une plasticité neuronale leur permettant, à l’issue d’une période de récupération, d’utiliser d’autres zones de leur cerveau pour effectuer la même tâche. Une telle plasticité est également observée dans le cas de patients privés d’un de leurs sens, les zones cérébrales habituellement affectées à ce sens pouvant alors être utilisées pour mieux utiliser les autres sens. L’étude de la plasticité neuronale à l’aide des potentiels évoqués pourrait permettre d’évaluer la progression du patient, ainsi que de mieux comprendre les mécanismes de cette plasticité, et l’efficacité d’un certain nombre de techniques de rééducation. Ce type d’application a été par exemple envisagé dans Neville et Bavelier [1996]; Posner et Abdullaev [1996]; Habib et al. [1996].

### 1.1.2 Activité électrique au niveau cellulaire.

Nous avons vu à la section précédente que les activités humaines telles que la perception, la mémorisation, le contrôle du fonctionnement de l’organisme (sommeil, sécrétions hormonales), font intervenir différentes zones cérébrales. Une manifestation de l’activité de ces zones est l’apparition d’un courant électrique membranaire dans les cellules nerveuses.
Tous les tissus vivants sont composés de cellules. Ces cellules sont spécialisées en fonction de leur rôle. N'importe quelle cellule peut présenter une différence de potentiel (le potentiel de repos) au niveau de sa membrane, mais la spécificité des cellules nerveuses et des cellules musculaires, lorsqu'elles sont excitées, est de pouvoir produire et propager des impulsions électriques (le potentiel d'action). Dans le cas des cellules musculaires, ce phénomène s'accompagne d'une contraction de la cellule. Dans d'autres cellules, ce phénomène peut également être associé à la fonction de la cellule.

Le cortex cérébral, ainsi que certains noyaux sous-corticaux, est constitué d'un grand nombre de cellules nerveuses, ou neurones, organisés selon différentes couches parallèles à la surface corticale. Ils sont composés d'un corps cellulaire, ou soma, se prolongeant par une longue fibre nerveuse, l'axone, et possédant de nombreuses ramifications, les dendrites. On trouve principalement dans le cortex deux types de neurones, les neurones pyramidaux et les neurones étoilés, qui diffèrent par leur forme (voir figure 1.1), mais aussi par leur organisation dans les couches corticales. Ainsi les neurones pyramidaux ont leurs dendrites orientés parallèlement entre elles et perpendiculairement à la surface corticale, tandis que les neurones étoilés n'ont pas d'orientation préférentielle.

La membrane d'une cellule nerveuse est formée d'une bicouche lipidique contenant des pores qui sous certaines conditions (incluant la différence de potentiel membranaire) peuvent être rendues perméables à certains ions qui leur sont spécifiques ($Na^+$, $K^+$, $Cl^-$). Des échanges d'ions ont lieu en permanence dans la cellule nerveuse au repos (via la "pompe à $Na^+/K^+$"), ce qui permet de maintenir l'équilibre à long terme de la composition ionique de la cellule, et donc le potentiel de repos. Des échanges peuvent également avoir lieu par les pores spécifiques, ou canaux ioniques, sous l'effet d'une excitation de la cellule, c'est ce qui engendre le potentiel d'action.

L'apparition, au delà d'un certain seuil, d'une différence de potentiel transmembranaire supérieure au potentiel de repos provoque l'ouverture de canaux ioniques qui vont rétablir une répartition homogène des ions de part et d'autre de la membrane. Les "pompes à $Na^+/K^+$" vont ensuite se remettre en action pour rétablir le potentiel de repos et la membrane va alors se surpolariser pour diminuer pendant un laps de temps l'excitabilité de la cellule (voir figure 1.2). Au cours de ce processus, la modification de l'équilibre ionique a de nouveau éloigné le potentiel transmembranaire du potentiel de repos, et l'onde électrique va ainsi se propager de proche en proche le long de l'axone. Une telle configuration de puits et de sources de courant peut être décrite, d'une manière très générale, par un développement multipolaire : la somme d'une composante monopolaire, d'une composante dipolaire, d'une composante quadripolaire, et ainsi de suite jusqu'à l'infini. Plus la dimension du multipôle est élevée et plus son influence sur le voisinage s'atténue rapidement avec la distance. Dans le cas de la propagation du potentiel d'action, l'entrée et la sortie à travers la membrane des ions $Na^+$ et
<figure><img src="image_9.png" /><figcaption>FIG. 1.1: Description d'un neurone pyramidal et d'un neurone étoilé. Une synapse excitatrice est représentée : elle provoque un puits de courant dans la région postsynaptique, ce qui provoque dans les milieux intra et extra-cellulaires des courants se refermant sur le reste de la membrane, lequel se présente donc comme une source de courant. Comme les sources de courant se situent préférentiellement vers le corps cellulaire, ces mouvements ioniques peuvent être modélisés en première approximation par un dipôle de courant, représenté par une flèche.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 1.2: Schéma de propagation du potentiel d'action</figcaption></figure>
$K^+$ peut être modélisée en première approximation par un quadripôle de courant, c'est à dire que les contributions de ses composantes monopolaire et dipolaire sont négligeables. Or les distributions de potentiel et de champ magnétique engendrées par un quadripôle s'annulent très vite avec la distance. *La contribution des potentiels d'action pour l'EEG et la MEG peut donc être négligée.*

Nous venons de voir que la propagation d'une onde électrique le long d'un axone se fait par simple conduction, mais le passage de l'information d'un neurone à un autre neurone fait intervenir un fonctionnement spécifique à travers les synapses. Il arrive parfois que l'espace synaptique soit suffisamment faible pour que le passage se fasse par simple conduction électrique, mais dans la plupart des cas, le passage de l'influx nerveux d'un neurone à l'autre fait appel à un processus biochimique particulier : lorsqu'une série de potentiels d'action parvient jusqu'au terminal présynaptique, cela provoque la libération de neuromédiateurs, lesquels en se fixant sur la membrane postsynaptique provoquent l'ouverture de canaux ioniques perméables aux ions $Na^+$ (synapse excitatrice) ou $Cl^-$ (synapse inhibitrice). L'entrée d'ions $Na^+$ en masse dans la cellule provoque l'apparition d'un potentiel d'action, tandis que l'entrée d'ions $Cl^-$ va abaisser la différence de potentiel transmembranaire et rendre plus difficile l'apparition d'un potentiel d'action au cas où des ions $Na^+$ viendraient à entrer également dans la cellule. Le soma du neurone est donc le lieu d'une sommation d'effets excitateurs et inhibiteurs permettant un traitement complexe de l'information véhiculée par les potentiels d'action. La distribution de puits et de sources de courant engendrée peut être modélisée en première approximation par un dipôle (voir figure 1.1). Cette composante du développement multipolaire subit une moins grande atténuation que le quadripôle de courant. Toutefois, il est communément admis [Nunez, 1981; Pernier et Bertrand, 1997] que *seuls les courants générés par l'activité synaptique synchrone d'une population de neurones de même orientation, comme les neurones pyramidaux perpendiculaires au cortex, en s'additionnant, engendrent un potentiel électrique et un champ magnétique mesurable à la surface de la tête, les autres sources de l'activité électrique cérébrale étant soumises à une atténuation telle qu'à la surface de la tête, leur effet peut être négligé.*

### 1.1.3 Modélisation de l'activité électrique cérébrale

Le problème qui nous concerne est de localiser les sources de l'activité électrique cérébrale à partir des données électroencéphalographiques et magnétoencéphalographiques. Il s'agit d'un problème inverse, dont la résolution passe par celle du problème direct associé : trouver les distributions de potentiel électrique et de champ magnétique qui seraient engendrées par des sources connues.
### 1.1.3.1 Modélisation des sources

Nous venons de voir que seule l'activité synaptique synchrone d'une population de neurones de même orientation pouvait engendrer une activité électrique mesurable à la surface de la tête. Nous avons également vu qu'en première approximation, l'activité synaptique pouvait être représentée par un dipôle de courant, et que de plus, les autres composantes du développement multipolaire engendrent une activité électrique qui s'atténue très rapidement avec la distance.

Quelle que soit l'étendue des sources, le problème peut se ramener à une distribution de sources dipolaires. Il suffit donc de savoir résoudre le problème direct avec ce type de modèle de source de base, le principe de superposition linéaire permettant de résoudre le cas général.

### 1.1.3.2 Modélisation de la géométrie de la tête

Le modèle de géométrie de la tête le plus simple<fnref n="1" /> et le plus largement utilisé est le modèle sphérique, [Rush et Driscoll, 1968] : la tête est assimilée à trois sphères homogènes de conductivités différentes (respectivement 0.45, 0.005625, et 0.45 Siemens/m), emboîtées les unes dans les autres et représentant le cerveau, l'os, et le scalp (voir figure 1.3). Dans ce cas, on peut calculer analytiquement [Ary et al., 1981; Sarvas, 1987] le potentiel et le champ magnétique engendrés par un dipôle donné. C'est pourquoi nous utiliserons ce modèle pour évaluer les performances de notre méthode de calcul.

L'une des activités de l'unité 280 de l'INSERM a consisté à développer et utiliser des modèles réalistes pour la géométrie de la tête. Il s'agit en pratique de modèles compartimentés comme dans le modèle sphérique, avec des conductivités homogènes par compartiment (certains comme l'os et le cerveau peuvent avoir une conductivité anisotrope), mais les frontières entre les compartiments sont obtenues à partir de piles d'IRM (Imagerie par Résonance Magnétique) fournissant une représentation anatomique d'une section de la tête. Un traitement des IRM permet d'extraire les contours du cerveau, de l'os, et du scalp. A partir de ces piles de contours, on génère des maillages de ces différentes surfaces (voir figure 1.4). C'est à ce modèle que devra s'appliquer notre méthode de calcul.

La résolution du problème direct avec ces modèles n'est réalisable qu'à l'aide de méthodes numériques. Différentes méthodes numériques peuvent être utilisées : la méthode des différences finies, la méthode des éléments finis, et la méthode des intégrales de frontière.

La résolution du problème direct par différences finies permet la prise en compte de l'anisotropie des différents milieux. La construction du maillage est facilitée car une grille régulière peut être utilisée. Elle se ramène à l'inversion d'une très grande matrice creuse. La

<footnote n="1">¹On trouve des modèles de géométrie encore plus simples comme le demi-espace [Cuffin et Cohen, 1977], mais ils ont plutôt été utilisés pour le magnétocardiogramme.</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 1.3: Exemple d'utilisation du modèle sphérique</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 1.4: Exemple de maillage réaliste de la tête</figcaption></figure>
définition d'un schéma aux différences finies pose quelques problèmes au niveau des surfaces de discontinuité, cependant la méthode a pu être mise en œuvre [Witwer et al., 1972; Marino et al., 1993; Hédou, 1997] de manière automatique à partir des IRM. L'inconvénient de cette méthode est qu'elle nécessite l'emploi d'ordinateurs massivement parallèles. Le temps de résolution d'un problème direct sur une Connexion Machine (partition CM5 de 32 nœuds) pour une grille cubique régulière de 256 éléments de coté était en 1997 de l'ordre de 10 mn d'après Hédou [1997].

On peut également prendre en compte l'anisotropie des milieux avec la méthode des éléments finis. Cette méthode permet de définir en chaque nœud du maillage tridimensionnel de la tête un tenseur de conductivité différent. Elle se ramène à l'inversion d'une grande matrice creuse. Appliquée pour la première fois à un modèle réaliste tridimensionnel de la tête par Thevenet et al. [1991]; Thevenet [1992]; Yan et al. [1991]; Marin [1997] et Marin et al. [1998], elle reste difficile à mettre en œuvre faute de procédure automatique permettant d'obtenir un maillage tridimensionnel de la tête incluant les surfaces de discontinuité. Le temps de résolution était en 1992 de l'ordre de 2 mn pour un maillage de 48000 tétraèdres avec un élément fini utilisant une interpolation linéaire sur des tétraèdres d'après Thevenet [1992].

La méthode que nous utiliserons est celle des intégrales de surface. Elle est également appelée méthode des éléments finis de surface, méthode des intégrales de frontière, ou méthode des éléments frontière. Elle se ramène à l'inversion d'une "petite" matrice pleine. Elle nécessite uniquement des maillages des différentes interfaces. Telle que nous l'appliquerons, elle ne nous permettra pas de prendre en compte les anisotropies. Cependant, il semblerait que celles-ci puissent être prises en compte à condition de modifier un peu l'approche [Zhou et Van Oosterom, 1994]. Les premiers à l'avoir appliquée à notre problème furent Hämäläinen [1987] et Meijs et al. [1987a]. Des algorithmes de maillage automatique des différentes surfaces ont pu être mis en œuvre [Fenollosa, 1993; Yvert, 1994; Mallet, 1996].

# 1.2 Formalisation mathématique du problème

## 1.2.1 Équations générales

Le potentiel et le champ magnétique engendrés par plusieurs dipôles de courant sont égaux à la somme des potentiels et des champs magnétiques créés par chacun de ces dipôles en vertu du principe de superposition. C'est pourquoi nous nous ramènerons à l'étude du potentiel et du champ magnétique engendrés par un seul dipôle.
### 1.2.1.1 Le potentiel électrique

L'équation de Poisson dans le milieu conducteur formé par la tête s'écrit :

$$ \nabla(\sigma\vec{\nabla}V) = \nabla(\vec{J}_p) \qquad (1.2.1.1) $$

où :

* $\vec{J}_p$ est la densité de courant due au dipôle de courant :

$$ \vec{J}_p = \vec{Q}\delta(\vec{r} - \vec{r}_o) \qquad (1.2.1.2) $$

où $\delta$ est le symbole de Dirac.

La densité de courant totale est la somme de ce cette densité de courant primaire avec les courants de conduction ohmique :

$$ \vec{J} = \vec{J}_p - \sigma \vec{\nabla} V \qquad (1.2.1.3) $$

* $\vec{Q}$ est le moment dipolaire du dipôle qui est placé en $\vec{r}_o$.

* $\sigma$ est la conductivité du milieu.

Par ailleurs, la continuité du potentiel et du courant ohmique sortant à travers les surfaces de discontinuité fixent les conditions aux limites. Notre problème peut donc s'écrire :

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

On a noté $[x]$ le saut de la variable $x$ en tout point d'une surface de discontinuité. On pourra également le noter $x^+ - x^-$.

### 1.2.1.2 Le champ magnétique

Le champ magnétique est donné par la loi de Biot et Savard :

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

Ce qui, combiné avec l'expression de la densité de courant (1.2.1.3), donne :

$$ \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' \qquad (1.2.1.6) $$

où $\vec{B}_o$ est le champ magnétique créé par le dipôle dans un milieu homogène infini :

$$ \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' \qquad (1.2.1.7) $$

Soit encore, avec l'équation (1.2.1.2),

$$ \vec{B}_o(\vec{r}) = \frac{\mu_o \vec{Q} \wedge (\vec{r} - \vec{r}_o)}{4\pi |\vec{r} - \vec{r}_o|^3} \qquad (1.2.1.8) $$
1.2.1.3 Modèle sphérique

Dans le cas d'un modèle sphérique à trois couches (scalp, os, cerveau), on peut calculer
analytiquement le potentiel [Ary et al., 1981] et le champ magnétique [Sarvas, 1987]. Les
équations, un peu compliquées, sont détaillées en annexe A.

Équation du potentiel électrique dans le modèle sphérique :

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

Équation du champ magnétique dans le modèle sphérique :

$$
\vec{B}(\vec{r}) = \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}
$$

1.2.2 Caractéristiques des champs de potentiel et des champs magnétiques

Les formules du champ magnétique dans le modèle sphérique impliquent plusieurs
conséquences qui sont essentielles pour comprendre la complémentarité des données MEG et
EEG.

Ainsi, dans le modèle sphérique en MEG :

* Un dipôle radial n'entraîne pas de champ magnétique.

* Le rayon de la sphère n'intervient pas.

* La décroissance du champ magnétique en fonction de la profondeur est importante, et un dipôle au centre de la sphère n'entraîne pas de champ magnétique.

* Les valeurs des conductivités n'interviennent pas.

**Remarque :** Dans le cas des modèle réalistes, la notion de dipôle radial ou tangential est moins bien définie. Fuchs et al. [1998a] ont proposé une définition des orientations quasi-radiale et quasi-tangentielles pour une position de dipôle donnée en se basant sur la décomposition en vecteurs propres du champ magnétique engendré (lead field) par les dipôles d'orientations parallèles aux axes du repère : soient $L_x$, $L_y$ et $L_z$ les champs magnétiques engendrés sur les $N_B$ voies MEG par les trois orientations selon les axes $x$, $y$ et $z$ du dipôle de moment unitaire. Pour un dipôle dont le moment a une orientation quelconque $\vec{M} = (M_x, M_y, M_z)$, le signal magnétique engendré est $B = L_x M_x + L_y M_y + L_z M_z$. Ce qui
peut s'écrire sous forme matricielle :

$$B = \begin{pmatrix} L_x & L_y & L_z \end{pmatrix} \begin{pmatrix} M_x \\ M_y \\ M_z \end{pmatrix} = LM$$

Le carré moyen du signal engendré sur les voies MEG est donc :

$$\frac{1}{N_B} B^t B = \frac{1}{N_B} M^t L^t LM$$

La matrice $L^tL$ est symétrique, donc diagonalisable, et le vecteur propre associé à la plus petite valeur propre donne l'orientation quasi-radiale du dipôle, tandis que les deux autres vecteurs propres correspondent aux orientations quasi-tangentielles. En général, la valeur propre correspondant à l'orientation quasi-radiale est très inférieure aux deux valeurs propres correspondant aux orientations quasi-tangentielles, qui elles, sont du même ordre de grandeur. Dans le modèle sphérique, l'orientation radiale du dipôle correspond au vecteur propre associé à la valeur propre zéro. On peut donc faire une analogie entre les orientations quasi-radiales dans le modèle réaliste, et les orientations radiales dans le modèle sphérique.

D'une manière générale, comme en première approximation la tête est globalement sphérique, la MEG enregistre moins bien que l'EEG les champs engendrés par les sources quasi-radiales et profondes. La contrepartie de cet inconvénient est que la MEG est fournit des topographies plus focales et en général l'influence des bruits physiologiques (activité musculaire, mouvements oculaires, clignements) sera plus locale. De même, elle est moins affectée que l'EEG par les valeurs des conductivités (voir annexe D).

D'autres aspects sont à prendre en compte pour comparer les avantages et les inconvénients des modalités MEG et EEG. L'installation des électrodes pour l'EEG est une opération d'autant plus longue que l'on utilise beaucoup d'électrodes. Il est plus rapide de placer un sujet dans sous un casque MEG, en revanche le sujet peut bouger légèrement par rapport aux voies MEG, tandis que les électrodes EEG, une fois posées, sont bien fixées sur le scalp. Par ailleurs, le coût d'une installation de MEG est très supérieur à celui d'une installation d'EEG.

### 1.2.3 La méthode des intégrales de surface

#### 1.2.3.1 Formulation intégrale

L'application de la deuxième formule de Green au problème (1.2.1.4) permet [Geselowitz, 1967] d'aboutir à l'équation intégrale de Fredholm (calculs détaillés en annexe B) :

$$(\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} \quad (1.2.1.9)$$
où

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

$V_o$ est le potentiel créé par le dipôle dans un milieu homogène infini de conductivité unitaire ; $\sigma_i^+$ et $\sigma_i^-$ sont les conductivités respectives à l'extérieur et à l'intérieur de la ième surface de discontinuité $S_i$ où se situe le point de calcul $r$ ; $\Omega_{i,r}^+$ et $\Omega_{i,r}^-$ sont les angles solides respectifs où le point $r$ voit l'extérieur et l'intérieur de la surface sur laquelle il se trouve.

Pour le champ magnétique, l'application du théorème de Gauss à l'équation (1.2.1.6) donne la formule de Geselowitz [1970] (calculs détaillés en annexe B) :

$$\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}' \quad (1.2.1.11)$$

### 1.2.3.2 Équations approchées matricielles

Le principe de la méthode des intégrales de surfaces peut être décrit comme suit (pour une information plus détaillée, voir Dautray et Lions [1985a] ou Brebbia et Dominguez [1989]).

L'équation que nous avons à résoudre peut se mettre sous la forme :

$$L(u) = b$$

où l'on cherche à connaître $u$ sur une surface $S$. L'opérateur $L$ est linéaire par rapport à $u$. Le principe de la méthode est de chercher à approcher $u$ par $\hat{u} \in W$, où $W$ est de dimension finie $n$ et engendré par une base de fonctions-test $\phi_k$ bien choisies :

$$\hat{u} = \sum_{k=1}^{n} \alpha_k \phi_k, \quad \alpha_k \in \mathbb{R}, \quad W = \left\{ \sum_{k=1}^{n} \beta_k \phi_k, \beta_k \in \mathbb{R} \right\}$$

Comme $\hat{u}$ n'est qu'une approximation de $u$, il subsiste un reste :

$$R(\hat{u}) = L(\hat{u}) - b \neq 0 \text{ à priori}$$

À partir de là, il existe principalement deux méthodes de résolution :

* Dans la méthode de Galerkin, on veut que le reste soit orthogonal à $W$. Ce qui donne les équations :

$$\forall k \in [1..n], \quad \int_S R(\hat{u}) \phi_k dS = 0.$$

qui équivaut à :

$$\forall k \in [1..n], \quad \sum_{j=1}^{n} \alpha_j \int_S L(\phi_j) \phi_k dS = \int_S b \phi_k dS$$
que l'on peut noter sous forme matricielle :

$$L_W \alpha = b_W,$$

où :

$$L_W = \left( \int_S L(\phi_k) \phi_j dS \right)_{1 \le j, k \le n}, \quad \alpha = (\alpha_k)_{1 \le k \le n}, \quad b_W = \left( \int_S b \phi_k dS \right)_{1 \le k \le n}$$

* Dans la méthode de collocation, on se donne $n$ points $\vec{x}_k$, $k \in [1..n]$ de la surface S, et on définit les fonctions-test de telle sorte que $\forall k, j \in [1..n]^2$, $\phi_k(\vec{x}_j) = \delta_{k,j}$. Il en résulte que $\forall k \in [1..n]$, $u(\vec{x}_k) = \alpha_k$. On veut alors que le reste soit tel que $\forall k \in [1..n]$, $R(\hat{u})(\vec{x}_k) = 0$. Ce qui donne les équations :

$$\forall k \in [1..n], \quad \sum_{j=1}^{n} \alpha_j L(\phi_j)(\vec{x}_k) = b(\vec{x}_k)$$

que l'on peut également noter sous forme matricielle :

$$L_W u = b_W,$$

où :

$$L_W = (L(\phi_k)(\vec{x}_j))_{1 \le j, k \le n}, \quad u = (u(\vec{x}_k))_{1 \le k \le n}, \quad b_W = (b(\vec{x}_k))_{1 \le k \le n}$$

Dans notre domaine, la méthode la plus utilisée pour résoudre notre problème est la méthode de collocation. C'est celle que nous avons utilisé dans ce travail.

Dans notre cas, nous appliquerons cette méthode aux équations de Fredholm (1.2.1.9) et de Geselowitz (1.2.1.11). Il nous faudra donc calculer des intégrales sur les surfaces de discontinuité. Pour cela, ces surfaces seront approchées par des polyèdres maillés. Deux ensembles de fonctions-test sont alors envisagés :

* dans l'approche centre de masse (COM) le potentiel est supposé constant sur chaque triangle du maillage :

$$\phi_k(\vec{x}) = \begin{cases} 1 & \text{si } \vec{x} \text{ est dans le triangle } k, \\ 0 & \text{sinon.} \end{cases}$$

Les nœuds $\vec{x}_k$ (points auxquels on calcule le potentiel) seront alors les centres de gravité des N triangles.

* dans l'approche interpolation linéaire (LINT) le potentiel varie linéairement entre les sommets des triangles du maillage :

$$\phi_k(\vec{x}) = \begin{cases} \frac{\det(\vec{x}, \vec{x}_l, \vec{x}_n)}{\det(\vec{x}_k, \vec{x}_l, \vec{x}_n)} & \text{si } \vec{x} \text{ est dans le triangle formé par } \vec{x}_k, \vec{x}_l, \vec{x}_n, \\ 0 & \text{sinon.} \end{cases}$$

Les nœuds $\vec{x}_k$ seront alors les sommets des triangles du maillage.
Ces approches permettent d'aboutir à des équations matricielles de la forme :

$$AV = V_o \qquad (1.2.1.12)$$

où les vecteurs $V$ et $V_o$ contiennent les valeurs de $V$ et $V_o$ aux nœuds du maillage pour l'EEG, et

$$B = B_o + MV \qquad (1.2.1.13)$$

où les vecteurs $B$ et $B_o$ contiennent les valeurs de $\vec{B}$ et $\vec{B}_o$ aux positions et selon l'orientation des spires pour la MEG. Le calcul des coefficients des matrices $A$ et $M$ est donné en annexe C.

### 1.2.3.3 La méthode de déflation [Lynn et Timlake, 1968]

La matrice $A$ du système linéaire obtenu par la méthode des intégrales de surface pour le potentiel électrique a une valeur propre nulle pour le vecteur propre constitué uniquement de $1 : \vec{e}^t = (11...11)$. Cela s'interprète d'un point de vue physique par le fait que le potentiel électrique est défini à une constante près. Le système $AV = V_o$ est donc un système singulier, et il faut y ajouter une condition supplémentaire pour rendre le système soluble.

Un artifice mathématique communément choisi est la déflation. Cette méthode consiste à remplacer la matrice $A$ par la matrice $C = A + \frac{\vec{e}^t \vec{e}}{N}$, où $N$ est la taille de la matrice, sous l'hypothèse que $0$ est valeur propre simple de $A$. On se retrouve avec le système inversible suivant : $CV = V_o$. La solution de ce système est alors solution du système $AV = V_o$, auquel on a ajouté la condition $\frac{V^t \vec{e}}{N} = 0$. La justification physique de cette condition est la suivante : on cherche le potentiel $V$ référencé de telle sorte que la somme des valeurs sur tous les nœuds des n couches soit nulle.

Dans le cas de l'approche centre de masse, cette méthode peut être encore perfectionnée car on connaît d'autres vecteurs propres de la matrice $A$ associés à des valeurs propres proches de 0. Cette méthode s'appelle la multidéflation. Elle reste applicable pour l'interpolation linéaire du potentiel [Yvert et al., 1996, p47].

### 1.2.3.4 L'approche problème isolé

Lors du développement de la méthode des intégrales de surface pour le problème de la localisation des sources de courant intracérébrales, Hämäläinen [1987] et Meijs et al. [1989] ont remarqué que la faible conductivité de l'os par rapport aux milieux voisins provoque d'importantes erreurs numériques. Ces erreurs numériques sont dues aux importants gradients de potentiels provoqués par la faible conductivité osseuse. Ils ont atténué ce problème en décomposant le potentiel $V$ en deux parties : la somme du potentiel $V_1$ que l'on obtient sur un
modèle à une couche (le cerveau) et du résidu $V_2$ lié aux autres couches : $V = V_1 + V_2$. La résolution se fait alors en deux étapes : la première consiste à calculer $V_1$ grâce à la résolution du système à une couche, et la deuxième étape consiste à calculer $V_2$ dans le modèle à 3 couches où l'on a modifié le second membre. Le potentiel $V_2$ a des gradients moins importants que $V$. Les équations sont les suivantes [Meijs et al., 1989] pour $\vec{r} \in \vec{S}_i$ :

$$ \sigma_3^- \Omega_{3,r}^- V_1(\vec{r}) = 4\pi V_o(\vec{r}) - \sigma_3^- \int_{S_3} V_1(\vec{r}') \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} \cdot d\vec{S} $$

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

$$ V_{o2}(\vec{r}) = -\delta_{i,3} \sigma_3^+ \Omega_{3,r}^+ V_1(\vec{r}) + \sigma_3^+ \int_{S_3} V_1(\vec{r}') \vec{\nabla} \frac{1}{|\vec{r} - \vec{r}'|} \cdot d\vec{S} $$

L'efficacité de cette méthode a été confirmée sur le problème direct par Meijs [Meijs et al., 1989] pour un dipôle dans un modèle sphérique à 4 couches, puis récemment sur le problème inverse par Zanow [Zanow, 1997, p27] pour dix dipôles tangentiels dans un modèle sphérique à 3 couches et des données de type MEG.

Le calcul des matrices pour l'approche problème isolé est détaillé en annexe E.

### 1.2.4 Méthodes de résolution du problème inverse

La résolution du problème inverse consiste à minimiser un critère mesurant l'écart entre le champ magnétique et/ou le potentiel électrique qui serait engendré(s) par une certaine configuration de sources de courant et ce qui a effectivement été mesuré. Dans le cas où un critère quadratique est choisi, il existe plusieurs méthodes pour trouver les paramètres qui minimisent ce critère.

#### 1.2.4.1 Choix d'un critère

##### Critère simple

Le choix le plus simple consiste à minimiser :

$$ S(\mathbf{d}) = \sum_{i=1}^{N_v} \sum_{t=1}^{N_{lat}} (B_{exp}^{i}(t) - B_{mod}^{i}(\mathbf{d}, t))^{2} $$

où d représente les paramètres (position et moment dipolaire) du (ou des) dipôle(s), $B_{exp}^{i}(t)$ le potentiel ou le champ magnétique mesuré à la voie i et à la latence t, et $B_{mod}^{i}(\mathbf{d}, t)$ le potentiel ou le champ magnétique que le dipôle de paramètres d engendre dans notre modèle à la voie i et à la latence t. $N_v$ est le nombre total de voies (MEG+EEG), et $N_{lat}$ est le nombre de latences.
On peut exprimer ce critère sous une forme matricielle :

$$S(\mathbf{d}) = \mathbf{f}^T \mathbf{f}.$$

où :

$$\mathbf{f} = \begin{pmatrix} B_{\text{exp}}^{1}(1) - B_{\text{mod}}^{1}(\mathbf{d}, 1) \\ B_{\text{exp}}^{1}(2) - B_{\text{mod}}^{1}(\mathbf{d}, 2) \\ \vdots \\ B_{\text{exp}}^{1}(N_{\text{lat}}) - B_{\text{mod}}^{1}(\mathbf{d}, N_{\text{lat}}) \\ B_{\text{exp}}^{2}(1) - B_{\text{mod}}^{2}(\mathbf{d}, 1) \\ \vdots \\ B_{\text{exp}}^{N_{\nu}}(N_{\text{lat}}) - B_{\text{mod}}^{N_{\nu}}(\mathbf{d}, N_{\text{lat}}) \end{pmatrix} = \begin{pmatrix} \mathbf{f}_{1,1} \\ \mathbf{f}_{1,2} \\ \vdots \\ \mathbf{f}_{1,N_{\text{lat}}} \\ \mathbf{f}_{2,1} \\ \vdots \\ \mathbf{f}_{N_{\nu},N_{\text{lat}}} \end{pmatrix}$$

### Critère pondéré

Le critère quadratique simple n'est pas très adapté lorsque, par exemple, le nombre de latences total est grand par rapport au nombre de latences où le dipôle a un moment dipolaire maximal. Pour toutes les latences où le moment dipolaire est faible, on va modéliser le bruit, ce qui n'est pas souhaitable. C'est pourquoi on peut pondérer le critère à minimiser de façon à ce que les latences de signal maximum aient plus d'importance dans le critère que les latences de signal minimum.

De même, dans le cas de données couplées MEG/EEG, on additionne des valeurs qui s'expriment dans des unités différentes. Il est préférable de normaliser les données à minimiser, par exemple en les pondérant par l'amplitude moyenne de chaque modalité.

Dans ce cas, on a toujours :

$$S(\mathbf{d}) = \mathbf{f}^T \mathbf{f}.$$

mais cette fois,

$$\mathbf{f}_{i,t} = w_{i,t} [B_{\text{exp}}^{i}(t) - B_{\text{mod}}^{i}(\mathbf{d}, t)]$$

où les $w_{i,t}$ représentent les poids de la pondération.

### Critère de maximum de vraisemblance

Le choix d'un critère pondéré dans le cas de données non homogènes met en évidence la difficulté à choisir le bon critère à minimiser. En effet le choix de ce critère a une grande influence sur le résultat de la minimisation. Intuitivement, on sent bien que le critère doit traiter de manière équitable toutes les informations élémentaires, sans qu'aucune ne soit privilégiée du fait de son amplitude plus grande. Or dans notre cas, ce qui est mesuré à une
électrode (ou à une spire) ne peut pas être très différent de ce qui est mesuré dans le voisinage. Notre information est corrélée spatialement. Cela signifie que l'information sur ce qui est mesuré à une électrode est en partie contenue également dans l'information mesurée au voisinage. Knuutila et Hämäläinen [1988], puis Lütkenhöner [1998] proposent de prendre en compte cette corrélation en utilisant un critère de maximum de vraisemblance.

Le principe est le suivant : pour le dipôle optimal d

$$B_{exp}^{i}(t) = B_{mod}^{i}(\mathbf{d}, t) + \delta_{i,t}$$

où $\delta$ est le bruit lié à la mesure et aux erreurs du modèle. Pour une autre réalisation du bruit, $B_{exp}$ aurait été différent. Supposons (c'est en fait une hypothèse très forte) que le bruit $\delta$ suive une loi gaussienne. La probabilité pour que $\delta$ se situe entre $(\delta_0 - \varepsilon)$ et $(\delta_0 + \varepsilon)$ s'écrit :

$$P = 2\varepsilon.e^{-\frac{1}{2}(\delta_0-\bar{\delta})^T \mathbf{W} (\delta_0-\bar{\delta})}$$

où $\mathbf{W}$ est l'inverse de la matrice de covariance de $\delta$, et $\bar{\delta}$ est la moyenne de $\delta$. Cette probabilité est maximale lorsque $(\delta_0 - \bar{\delta})^T \mathbf{W} (\delta_0 - \bar{\delta})$ est minimal, ce qui permet de proposer comme critère à minimiser :

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

où $\mathbf{W}$ est une matrice symétrique définie positive qui représente l'inverse d'un estimateur de la matrice de covariance, et :

$$\mathbf{f}_{i,t} = [B_{exp}^{i}(t) - B_{mod}^{i}(\mathbf{d}, t)] - \overline{[B_{exp}^{i}(t) - B_{mod}^{i}(\mathbf{d}, t)]}$$

Notons que dans le cas où la moyenne du bruit est nulle, et où $\mathbf{W}$ est diagonale (c'est-à-dire que le bruit est indépendant sur chaque voie et chaque latence), on se ramène au cas du critère pondéré, et que dans le cas où la variance du bruit est identique pour chaque voie et chaque latence, $\mathbf{W} = \frac{1}{\sigma^2} \mathbf{I}$, on se ramène au cas du critère quadratique simple.

De plus, comme $\mathbf{W}$ est symétrique définie positive, $\exists \mathbf{R}, \mathbf{W} = \mathbf{R}^T \mathbf{R}$. En posant $\mathbf{g} = \mathbf{R} \mathbf{f}$, on a alors :

$$S(\mathbf{d}) = \mathbf{g}^T \mathbf{g},$$

ce qui nous ramène au critère simple, à une modification linéaire des données près.

### 1.2.4.2 Méthodes de minimisation

#### La méthode de plus profonde descente [Cauchy, 1847]

On souhaite minimiser $S(\mathbf{d}) = \mathbf{f}^T \mathbf{W} \mathbf{f}$. Les paramètres des dipôles sont initialisés à $\mathbf{d}_0$. Ensuite, à chaque itération, connaissant $\mathbf{d}_i$, on calcule le gradient

$$\frac{\partial S}{\partial \mathbf{d}}(\mathbf{d}_i) = 2\mathbf{v},$$
et on cherche $t$ tel que

$$S(\mathbf{d}_i + t\mathbf{v}) < S(\mathbf{d}_i).$$

Si un tel $t$ ne peut être trouvé, c'est que l'on se situe dans un minimum local. Sinon, on prend $\mathbf{d}_{i+1} = \mathbf{d}_i + t\mathbf{v}$, et on recommence les itérations. On peut montrer que cette méthode converge lorsque la fonction $S(\mathbf{d})$ est convexe. Lorsque ce n'est pas le cas, la méthode converge vers un minimum local ou un point-selle. Dans tous les cas, la convergence peut être très lente.

## La méthode de Gauss-Newton

On souhaite trouver $\mathbf{d}$ tel que $\mathbf{v}(\mathbf{d}) = 0$. Les paramètres des dipôles sont initialisés à $\mathbf{d}_0$. Ensuite, à chaque itération, connaissant $\mathbf{d}_i$, on calcule $\mathbf{v}(\mathbf{d}_i)$. On peut écrire :

$$2\mathbf{v}_j(\mathbf{d}_i) = \frac{\partial \mathbf{f}^T}{\partial d_j} \mathbf{W} \mathbf{f} + \mathbf{f}^T \mathbf{W} \frac{\partial \mathbf{f}}{\partial d_j}$$

et comme $\mathbf{W}$ est symétrique,

$$\mathbf{v}(\mathbf{d}_i) = \mathbf{J}^T \mathbf{W} \mathbf{f}$$

où :

$$\mathbf{J}_{k,j} = \frac{\partial \mathbf{f}_k}{\partial d_j}.$$

Par ailleurs,

$$\mathbf{v}_j(\mathbf{d}_i + \mathbf{q}) = \mathbf{v}_j(\mathbf{d}_i) + \sum_{k=1}^{n} q_k \frac{\partial v_j(\mathbf{d}_i)}{\partial d_k} + O(q^2).$$

Or

$$\frac{\partial v_j(\mathbf{d}_i)}{\partial d_k} = (\mathbf{J}^T \mathbf{W} \mathbf{J})_{j,k} + (\frac{\partial \mathbf{J}^T}{\partial d_k})_{j,k} \mathbf{W} \mathbf{f},$$

d'où :

$$\mathbf{v}(\mathbf{d}_i + \mathbf{q}) = \mathbf{J}^T \mathbf{W} \mathbf{f} + \mathbf{J}^T \mathbf{W} \mathbf{J} \mathbf{q} + \text{termes en } \frac{\partial^2 \mathbf{f}}{\partial d_j \partial d_k} + O(q^2).$$

On néglige les termes en $q^2$ et les dérivées secondes (ce qui est rarement justifié). La solution est alors donnée quand

$$\begin{aligned} 0 &= \mathbf{J}^T \mathbf{W} \mathbf{f} + \mathbf{J}^T \mathbf{W} \mathbf{J} \mathbf{q} \\ \mathbf{q} &= -(\mathbf{J}^T \mathbf{W} \mathbf{J})^{-1} \mathbf{J}^T \mathbf{W} \mathbf{f} \end{aligned}$$

En pratique, $\mathbf{d}_{i+1} = \mathbf{d}_i + \mathbf{q}$ n'est pas exactement la solution, mais on réitère le procédé jusqu'à ce que $\mathbf{q}$ soit inférieur à une valeur minimale prédéfinie ou que $S(\mathbf{d}_i + \mathbf{q}) \ge S(\mathbf{d}_i)$.

Cette méthode échoue souvent car les hypothèses sur les dérivées secondes sont trop fortes. De plus, la matrice $\mathbf{J}^T \mathbf{W} \mathbf{J}$ peut poser des problèmes de régularité.
## Méthode de Marquardt [1963]

La méthode de Marquardt est une méthode intermédiaire entre la méthode de plus profonde descente et la méthode de Gauss-Newton. Il propose de chercher les solutions de :

$$ (\mathbf{J}^T \mathbf{W} \mathbf{J} + \lambda \mathbf{I}) \mathbf{q} = -\mathbf{J}^T \mathbf{W} \mathbf{f} \quad (1.2.1.14) $$

Comme $v = \mathbf{J}^T \mathbf{W} \mathbf{f}$, si $\lambda$ est très grand par rapport à la norme de $\mathbf{J}^T \mathbf{W} \mathbf{J}$, on se trouve dans le cas d'une méthode de plus profonde descente, et si $\lambda = 0$, on retrouve la méthode de Gauss-Newton.

Afin de s'affranchir également de problèmes d'échelle entre les paramètres, il propose d'utiliser à la place de $d$,

$$ d' = Dd $$

où $D$ est une matrice diagonale inversible ayant tous ses coefficients strictement positifs. La matrice jacobienne contenant les dérivées de $f$ par rapport aux $d'_k$ devient alors :

$$ \mathbf{J}' = \mathbf{J} \mathbf{D}^{-1} $$

L'équation 1.2.1.14 s'écrit alors :

$$ \begin{align} ((\mathbf{J}')^T \mathbf{W} \mathbf{J}' + \lambda \mathbf{I}) \mathbf{q}' &= -(\mathbf{J}')^T \mathbf{W} \mathbf{f} \nonumber \\ (\mathbf{D}^{-1} \mathbf{J}^T \mathbf{W} \mathbf{J} \mathbf{D}^{-1} + \lambda \mathbf{I}) \mathbf{D} \mathbf{q} &= -\mathbf{D}^{-1} \mathbf{J}^T \mathbf{W} \mathbf{f} \nonumber \\ (\mathbf{J}^T \mathbf{W} \mathbf{J} + \lambda \mathbf{D}^2) \mathbf{q} &= -\mathbf{J}^T \mathbf{W} \mathbf{f} \tag{1.2.1.15} \end{align} $$

Marquardt propose de suivre la suggestion de Levenberg [1944] :

$$ \mathbf{D}_{ii} = (\mathbf{J}^T \mathbf{W} \mathbf{J})_{ii} $$

Nash [1979] suggère plutôt

$$ \mathbf{D}_{ii} = (\mathbf{J}^T \mathbf{W} \mathbf{J})_{ii} + 1, $$

ce qui permet d'éviter que l'un des $\mathbf{D}_{ii}$ soit nul ou même numériquement trop proche de 0. Dans ce cas, la matrice $(\mathbf{J}^T \mathbf{W} \mathbf{J} + \lambda \mathbf{D}^2)$ est symétrique définie positive pourvu que $\lambda$ soit assez grand, et on peut donc résoudre l'équation 1.2.1.15 par décomposition de Choleski.

Marquardt [1963] propose de choisir $\lambda = 0.1$ au départ, puis tant que $S(\mathbf{d}_i + \mathbf{q}) < S(\mathbf{d}_i)$, multiplier $\lambda$ par 0.1, et passer à l'itération suivante avec $\mathbf{d}_{i+1} = \mathbf{d}_i + \mathbf{q}$; en revanche, si $S(\mathbf{d}_i + \mathbf{q}) \ge S(\mathbf{d}_i)$, multiplier $\lambda$ par 10, et passer à l'itération suivante en gardant le même $\mathbf{d}_i$. Nash [1979] propose de choisir plutôt comme facteurs 0.4 et 10, et d'initialiser à 0.0001, mais la méthode reste la même.
<figure><img src="image_13.png" /><figcaption>FIG. 1.5: Algorithme de Marquardt [1963] modifié par Nash [1979]</figcaption></figure>

### Prise en compte de la dépendance linéaire au moment dipolaire

Le champ magnétique et le potentiel électrique engendrés par un dipôle de courant dépendent linéairement de l'amplitude du moment dipolaire. Or dans le cas d'une dépendance linéaire, on peut calculer directement, sans passer par un processus itératif de minimisation, les valeurs des paramètres qui optimisent un critère quadratique. On pourra alors améliorer l'algorithme de Marquardt en séparant les paramètres correspondant à l'amplitude du moment dipolaire des paramètres de position [Scherg, 1990]. Seuls les paramètres de position et d'orientation restent libres pour l'algorithme de Marquardt, les moments dipolaires sont déduits à chaque itération par optimisation linéaire.

### Condition nécessaire au minimum

Notons $d_i$ les paramètres de position et d'orientation des dipôles, et $M_{i,t}$ les moments dipolaires pour chaque dipôle $i$ et pour chaque latence $t$, et $B^e(d_i)$ le potentiel ou le champ magnétique modèle à la voie $e$ pour un dipôle de moment dipolaire égal à 1. On a alors :

$$B_{\text{mod}}^{e}(\mathbf{d}, t) = \sum_{i=1}^{N_{\text{dip}}} B^{e}(d_{i})M_{i,t}$$

Ce que l'on peut noter sous la forme matricielle :

$$B_{\text{mod}}(\mathbf{d}) = B(\mathbf{d})M$$
où

$$
B_{\text{mod}}(\mathbf{d}) = \left( \begin{array}{c} B^1(\mathbf{d}, t_1) \\ \vdots \\ B^1(\mathbf{d}, t_{N_{\text{lat}}}) \\ B^2(\mathbf{d}, t_1) \\ \vdots \\ B^{N_{\text{voies}}}(\mathbf{d}, t_{N_{\text{lat}}}) \end{array} \right) \Bigg|_{N_{\text{voies}}*N_{\text{lat}}} B(\mathbf{d}) = \begin{array}{c} \left( \begin{array}{cccc} \tilde{B}(\mathbf{d}) & 0 & \dots & 0 \\ 0 & \tilde{B}(\mathbf{d}) & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \tilde{B}(\mathbf{d}) \end{array} \right) \end{array} \Bigg|_{N_{\text{voies}}*N_{\text{lat}}}
$$

$$
\tilde{B}(\mathbf{d}) = \left( \begin{array}{ccc} B^1(\mathbf{d}_1) & \dots & B^1(\mathbf{d}_{N_{\text{dip}}}) \\ \vdots & \ddots & \vdots \\ B^{N_{\text{voies}}}(\mathbf{d}_1) & \dots & B^{N_{\text{voies}}}(\mathbf{d}_{N_{\text{dip}}}) \end{array} \right) \Bigg|_{N_{\text{voies}}} \quad \mathbf{M} = \left( \begin{array}{c} M_{1,t_1} \\ \vdots \\ M_{1,t_{N_{\text{lat}}}} \\ M_{2,t_1} \\ \vdots \\ M_{N_{\text{dip}},t_{N_{\text{lat}}}} \end{array} \right) \Bigg|_{N_{\text{dip}}*N_{\text{lat}}}
$$

On souhaite minimiser, les paramètres de position d étant fixes,

$$
\begin{align*}
S(\mathbf{M}) &= (B_{\text{exp}} - B_{\text{mod}}(\mathbf{d}))^T W (B_{\text{exp}} - B_{\text{mod}}(\mathbf{d})) \\
&= (B_{\text{exp}} - B(\mathbf{d})\mathbf{M})^T W (B_{\text{exp}} - B(\mathbf{d})\mathbf{M})
\end{align*}
$$

On a alors, comme W est symétrique,

$$
\nabla(S(\mathbf{M})) = 2 \mathbf{B}(\mathbf{d})^T W (\mathbf{B}_{\text{exp}} - \mathbf{B}(\mathbf{d})\mathbf{M})
$$

Lorsque le minimum est atteint, le Jacobien est nul, on a donc :

$$
\mathbf{B}(\mathbf{d})^T W (\mathbf{B}_{\text{exp}} - \mathbf{B}(\mathbf{d})\mathbf{M}) = 0
$$

donc

$$
\mathbf{B}(\mathbf{d})^T W \mathbf{B}(\mathbf{d}) \mathbf{M} = \mathbf{B}(\mathbf{d})^T W \mathbf{B}_{\text{exp}} \quad (1.2.1.16)
$$

Supposons maintenant que M annule le Jacobien. On a alors :

$$
\begin{align*}
\forall \partial \mathbf{M}, S(\mathbf{M} + \partial \mathbf{M}) &= (B_{\text{exp}} - B(\mathbf{d})(\mathbf{M} + \partial \mathbf{M}))^T W (B_{\text{exp}} - B(\mathbf{d})(\mathbf{M} + \partial \mathbf{M})) \\
&= (B_{\text{exp}} - B(\mathbf{d})\mathbf{M})^T W (B_{\text{exp}} - B(\mathbf{d})\mathbf{M}) \\
&\quad -2 \partial \mathbf{M}^T B(\mathbf{d})^T W (B_{\text{exp}} - B(\mathbf{d})\mathbf{M}) \\
&\quad + (B(\mathbf{d})\partial \mathbf{M})^T W (B(\mathbf{d})\partial \mathbf{M}) \\
&= S(\mathbf{M}) + 0 + (B(\mathbf{d})\partial \mathbf{M})^T W (B(\mathbf{d})\partial \mathbf{M}) \\
&\geq S(\mathbf{M}) \text{ car } W \text{ est positive}
\end{align*}
$$
Donc si M est solution de 1.2.1.16, alors c'est un minimum.

### Existence et unicité du minimum

Si $B(d)^T W B(d)$ est inversible, il existe une unique solution à 1.2.1.16. Supposons que ce
ne soit pas le cas, alors $\exists M \neq 0$, $B(d)^T W B(d) M = 0$, et donc $M^T B(d)^T W B(d) M = 0$, et
comme W est définie positive, $B(d)M = 0$. En pratique, ce cas peut arriver, par exemple dans
le modèle sphérique en MEG, lorsqu'un dipôle est radial ou au centre de la sphère, ou lorsque
deux dipôles ont la même position et des moments dipolaires opposés. Ces cas doivent être
evités en se fixant des contraintes sur la position et l'orientation du dipôle.

Soit $W = R^T R$ la décomposition de Choleski de W. L'équation 1.2.1.16 peut s'écrire :

$$B(d)^T R^T R B(d) M = B(d)^T R^T R B_{exp}$$

Ce qui équivaut, en notant $Z = RB(d)$ et $Y = RB_{exp}$, à :

$$Z^T Z M = Z^T Y$$

Or cette classe de problèmes admet toujours au moins une solution (voir la démonstration page 188 de l'annexe H). Toute solution composée de cette solution plus un élément du noyau de $B(d)$ est également solution.

En résumé, il existe au moins une solution à 1.2.1.16, mais elle n'est unique que si $B(d)^T W B(d)$ est inversible, ce qui équivaut à dire qu'elle est symétrique définie positive.
Dans tous les cas, on peut utiliser la décomposition de Choleski de $B(d)^T W B(d)$ en prenant des précautions en cas de singularité pour trouver une solution à 1.2.1.16 orthogonale au noyau de $B(d)$ (voir annexe H).

### 1.2.5 Mise en œuvre

Dans notre laboratoire, la solution logicielle mise en place pour résoudre le problème inverse est la suivante :

* une bibliothèque de fonctions "BIBMODEL" permet le calcul du problème direct. Avant mon arrivée au laboratoire, seuls les problèmes direct EEG sphérique et réalistes par la méthode des intégrales de surface avec l'approche COM étaient programmés. J'y ai ajouté en EEG l'approche LINT et en MEG la résolution du problème direct avec les modèles sphériques et réalistes.

* Le logiciel "ERPA" est une solution complète permettant, au travers d'une interface homme-machine intuitive, de manipuler des fichiers de potentiels évoqués afin de tracer des courbes ou des cartes d'une part, et afin de préparer la résolution du problème inverse d'autre part en choisissant le type de modèle désiré, le nombre de sources, les
conductivités des milieux. J'y ai ajouté la possibilité de lire des fichiers de type MEG contenant les informations nécessaire à la résolution du problème inverse : en particulier les positions des spires et la description des combinaisons linéaires de mesures par spires permettant d'aboutir aux mesures par voie. J'y ai également ajouté l'interface nécessaire pour pouvoir choisir l'approche LINT, la résolution du problème inverse en MEG sphérique ou réaliste, ou par couplage MEG/EEG.

* Le logiciel "PEMODEL" est dédié à la résolution du problème inverse. C'est un logiciel interactif qui permet de visualiser en temps réel l'évolution de la résolution du problème : nombre d'itération, position et orientation des sources, valeur du critère. Une option permet également la résolution du problème inverse de manière non-interactive ("batch") mais plus rapide.

Pour résoudre le problème inverse, "PEMODEL" offre la possibilité de mettre des bornes et des contraintes sur les paramètres. Cela permet d'une part de restreindre l'espace de recherche, et d'autre part d'empêcher les dipôles de prendre des valeurs aberrantes (par exemple en dehors de la tête). Certaines contraintes diminuent le nombre de degrés de libertés pour les paramètres : par exemple, les contraintes de symétrie entre les sources, la contrainte de tangentialité pour l'orientation. La minimisation est alors faite en faisant varier les paramètres libres, et en calculant les paramètres liés à chaque itération. Les autres contraintes fixent des bornes pour les paramètres libres. Lorsque le dipôle se trouve sur une borne, tous les déplacements ne peuvent pas être autorisés : dans le cas où le signe de sa dérivée amènerait un paramètre libre à se déplacer hors de ses bornes, on le bloque; dans le cas contraire on le débloque, ce qui permet au dipôle de ne pas rester coincé lorsqu'une borne est atteinte en cours de minimisation.

Le programme "PEMODEL" doit permettre de choisir entre les critères simples, pondérés, et de maximum de vraisemblance. En choisissant de se donner une matrice de poids W, cet objectif peut facilement être atteint en utilisant la matrice identité, une matrice diagonale ou une estimation de l'inverse de la matrice de covariance. Cette matrice peut soit être donnée par le chercheur, soit être calculée par le programme. Le programme "PEMODEL" distingue les cas où la matrice de poids est diagonale et où elle est pleine de façon à accélérer les calculs.

J'ai ajouté à "PEMODEL" les appels nécessaires lorsque l'approche LINT est utilisée, ainsi que la résolution en MEG sphérique et réaliste et le couplage MEG/EEG, et la prise en compte de critères définis par une matrice de poids.

* Le logiciel "VISIMAGE" permet la visualisation d'IRM anatomiques. Il peut lire les fichiers de résultats de "PEMODEL" pour superposer les solutions trouvées avec les images anatomiques. En amont du problème inverse, il permet le choix d'une sphère
lorsque ce modèle a été choisi, ou de visualiser le résultat de modules de segmentation semi-automatiques pour préparer la construction des maillages de surface lorsque la modélisation réaliste a été choisie. Ce logiciel est en cours de migration vers un logiciel plus vaste, “ACTIVIS”. Je n’ai pas eu besoin d’y ajouter de modification dans le cadre de cette thèse.

* Le logiciel "MARS" permet la construction semi-automatique des maillages de surface à partir de piles de contours obtenues par segmentation des IRM. C'est un logiciel interactif à l'interface intuitive. La création d'un arbre de connection entre les contours des différentes surfaces est automatique, toutefois en cas d'erreur de l'euristique on peut corriger manuellement à travers une interface graphique cet arbre de connection. A partir de cet arbre, le maillage est réalisé automatiquement à partir du pas de discrétisation défini sur chaque contour. Les intersections éventuelles sont détectées, et une correction manuelle peut être effectuée à travers l'interface graphique. Ensuite le maillage peut être affiné localement ou bien simplifié. Ce programme permet également une cartographie réaliste des potentiels ou des champs magnétiques. Je n'ai pas eu besoin d'y ajouter de modification dans le cadre de cette thèse.

* Le programme "BEM" est une interface avec la bibliothèque "BIBMODEL" qui permet de calculer les matrices et leurs inverses ou leurs décompositions LU pour la méthode des intégrales de surface. Il permet également la génération de fichiers de potentiels et de champ magnétique simulés, obtenus avec un modèle réaliste ou sphérique. J'y ai ajouté la possibilité de choisir l'approche LINT, et le modèle sphérique ou réaliste en MEG.
## 1.3 Synthèse bibliographique de l'évaluation des modèles

Cette section présente une revue bibliographique des évaluations des performances de différentes approches permettant de résoudre le problème inverse et datant d'avant le début de cette thèse. Ces études sur les effets du choix d'un modèle réaliste ou non, sur l'apport de l'approche problème isolé, l'apport de l'affinage local des maillages, l'effet du choix de la fonction d'interpolation pour les modèles réalistes, l'apport du couplage des données MEG et EEG permettront de dresser un plan de travail pour cette thèse en fonction de l'état de l'art.

### 1.3.1 Apport des modèles réalistes par rapport aux modèles sphériques

#### 1.3.1.1 Apport des modèles réalistes pour l'EEG

**Introduction des modèles réalistes**

Rush et Driscoll [1968] ont introduit le modèle sphérique à trois couches permettant de mieux rendre compte qu'un modèle sphérique homogène de la faible conductivité de l'os. L'influence des paramètres de ce modèle (rayons des sphères, nombre de couches, conductivités des différentes couches) a alors été étudiée.

Ary et al. [1981] ont placé 15 dipôles radiaux et 15 dipôles tangentiels dans des modèles sphériques à trois couches avec 5 valeurs différentes pour l'épaisseur du scalp, 2 valeurs différentes pour celle de l'os et 3 valeurs différentes pour le rapport entre la conductivité de l'os et celle du scalp et du cerveau. Ils ont ensuite résolu le problème inverse EEG avec un modèle sphérique homogène. L'erreur sur l'excentricité du dipôle était d'autant plus importante que l'excentricité et les épaisseurs du scalp et de l'os étaient importantes. Ils ont ainsi montré qu'une erreur raisonnable sur les épaisseurs des couches pouvait introduire une erreur de 7% sur l'estimation de l'excentricité tandis qu'une erreur sur la conductivité de l'os n'introduisait qu'une erreur de 2%. Cette étude montre l'importance de la qualité de la description géométrique du modèle pour la résolution du problème inverse. C'est dans ce contexte que des modèles réalistes de la tête ont été introduits [Hämäläinen, 1987; Meijs et al., 1987a; 1989].

**Évaluation des modèles réalistes par simulation**

Cuffin [1990] a placé trois dipôles (1 position et 3 orientations) dans 8 modèles différents de la tête : un modèle sphérique à trois couches pour lequel la résolution du problème direct est analytique, et sept modèles de type intégrale de surface (BEM) à trois couches avec des maillages de surface composés de 1088 triangles par couche, l'un sphérique, cinq ellipsoïdaux et le dernier sphérique pour les deux couches externes et demi-sphérique pour la couche interne. Il ne s'agit donc pas tout à fait de l'évaluation des modèles réalistes, mais plutôt de
l'évaluation de l'effet de la forme du modèle. Le problème direct a été résolu en EEG et en MEG pour chaque modèle, et le problème inverse a alors été résolu à partir de ces données simulées avec un modèle sphérique à trois couches. Pour le modèle sphérique analytique en EEG, le dipôle était parfaitement retrouvé. Avec le modèle BEM sphérique, l'erreur de localisation en EEG restait inférieure à 2 mm. Avec les modèles BEM ellipsoïdaux, l'erreur de localisation en EEG variait entre 2.3 mm et 17.7 mm. Avec le modèle BEM où le cerveau était modélisé par une demi-sphère, l'erreur de localisation en EEG variait entre 2.9 mm et 4.3 mm.

Roth et al. [1993] ont placé 24 dipôles (8 positions et 3 orientations axiales par position) dans un modèle réaliste de la tête composé de trois couches avec 1600 triangles par couche. Le potentiel électrique pour 21 positions d'électrodes a ainsi été calculé avec la méthode des intégrales de surface. Ensuite, le problème inverse a été résolu avec un modèle sphérique à trois couches et avec des modèles réalistes différents de celui utilisé pour le problème direct et composés de 256, 576, 1024 et 1600 triangles par couche. Les erreurs de localisation avec le modèle sphérique étaient en moyenne de 2 cm alors que les modèles réalistes conduisaient à des erreurs de localisation moyennes de 0.8 cm avec 256 ou 576 triangles par couche et 0.4 cm avec 1024 ou 1600 triangles par couche.

Zanow et Peters [1995] ont placé deux dipôles à deux profondeurs différentes de 5 mm et 23 mm dans un modèle réaliste à trois couches avec environ 800 triangles par couches. Le potentiel électrique pour 52 positions d'électrodes a ainsi été calculé avec la méthode des intégrales de surface. Le problème inverse a alors été résolu avec trois modèles sphériques différents à trois couches : l'un ajusté au scalp, l'autre aux électrodes et le dernier au cerveau. Deux modèles réalistes ont également été utilisés pour résoudre le problème inverse, l'un était celui utilisé pour le problème direct et l'autre était différent mais avait la même densité de triangles par couche (il avait été obtenu en recommençant la procédure d'obtention du maillage à partir des IRM). Le problème inverse a été résolu en ajoutant 0%, 5% et 10% de bruit blanc gaussien. Le modèle sphérique qui donnait les meilleurs résultats était celui ajusté aux électrodes. L'erreur de localisation avec ce modèle était de 4.5 mm et 6.5 mm en l'absence de bruit, 5.1 mm et 7.3 mm en moyenne avec 5% de bruit et 5.7 mm et 10.6 mm en moyenne avec 10% de bruit. Les erreurs de localisation étaient d'autant plus importantes que le niveau de bruit et la profondeur du dipôle augmentaient. Le modèle réaliste différent de celui du problème direct donnait de meilleurs résultats que le modèle sphérique en l'absence de bruit (1.3 mm et 4 mm) et lorsque le niveau de bruit était de 5% (2.2 mm et 6.6 mm) mais quand le niveau de bruit était de 10% l'erreur de localisation était plus importante pour le dipôle profond (3.6 mm et 10.9 mm).

Yvert et al. [1997] ont réalisé une simulation très systématique en plaçant 2000 dipôles sur la surface du cortex et orientés perpendiculairement à celle-ci dans un modèle réaliste à trois couches composé de 4000 triangles environ et localement affiné dans le cas des dipôles
périphériques. Pour chaque dipôle le potentiel électrique était calculé en 19, 32 et 63 positions d'électrodes. Le problème inverse a alors été résolu avec différents modèles sphériques et un modèle réaliste différent de celui utilisé pour le problème direct. Ils ont trouvé une erreur de localisation moyenne avec le meilleur modèle sphérique de 5-6 mm dans la partie supérieure de la tête et 15-25 mm dans la partie inférieure de la tête. L'utilisation d'un modèle réaliste uniformément maillé permet d'aboutir à une erreur de localisation moyenne de 4.4 mm.

## Évaluation des modèles réalistes avec des dipôles implantés chez des patients épileptiques

Dans deux articles traitant les mêmes données EEG [Cuffin et al., 1991; Cuffin, 1996], 28 dipôles implantés chez trois patients épileptiques ont été localisés en utilisant d'une part un modèle sphérique à quatre couches (le scalp, l'os, le liquide céphalo-rachidien et le cerveau) et d'autre part un modèle réaliste composé de trois couches avec 1088 triangles par couche. Le modèle sphérique conduisait à des erreurs de localisation de 11 mm en moyenne tandis que les modèles réalistes permettaient d'améliorer cette erreur de manière significative chez un des trois patients. Lorsque le rapport signal sur bruit était faible, chez les deux autres patients, le modèle réaliste n'apportait pas toujours d'amélioration de la localisation.

Homma et al. [1994] ont recueilli 19 enregistrements de dipôles implantés chez des patients épileptiques et 5 enregistrements avec deux dipôles. Ils ont utilisé un modèle réaliste de la tête constitué de trois couches avec 2496, 1104 et 2208 triangles respectivement pour le scalp, l'os et le cerveau, et ont résolu le problème inverse avec la méthode des intégrales de surface et une interpolation linéaire du potentiel. Lorsqu'un seul dipôle était présent, ils ont trouvé une erreur de localisation moyenne de 8.5 mm tandis que lorsque deux dipôles étaient présents, l'erreur moyenne était de 18 mm. Le modèle sphérique n'a pas été utilisé dans cet article.

### 1.3.1.2 Apport des modèles réalistes pour la MEG

Il est fréquent de lire que la MEG est moins affectée que l'EEG par la forme exacte de la géométrie de la tête et par les valeurs des conductivités. Un argument souvent évoqué pour justifier l'intérêt de la MEG par rapport à l'EEG est que la faible conductivité de l'os rend l'EEG moins focale et donc moins apte à détecter des sources occupant de petites régions<fnref n="2" /> [Cuffin et Cohen, 1979]. On lit également que comme la tête est en première approximation sphérique, la MEG est globalement moins sensible que l'EEG aux variations de conductivité

<footnote n="2">²Extrait de Cuffin et Cohen [1979] : "These results suggest that the MEG may detect a more restricted group of sources than the EEG, i.e., sources which occupy a smaller region..."</footnote>
[Cuffin et Cohen, 1979] ou que l'effet des courants de conduction est négligeable en MEG [Hämäläinen et Sarvas, 1989].

### Évaluation des modèles réalistes par simulation

Dans la simulation de Cuffin [1990] décrite à la section précédente, pour trois dipôles (1 position et 3 orientations), le problème direct réaliste était résolu en MEG dans 7 modèles différents de la tête : un modèle sphérique BEM, cinq ellipsoïdaux BEM et le dernier sphérique pour le scalp et l'os et demi-sphérique pour le cerveau. Le problème inverse était résolu avec un modèle sphérique analytique. L'erreur de localisation en MEG avec le modèle sphérique BEM était inférieure à 0.1 mm pour les dipôles tangentiels, et un dipôle de moment nul était retrouvé pour le dipôle radial. Avec les modèles BEM ellipsoïdaux et les dipôles tangentiels l'erreur de localisation variait entre 2.6 mm et 12.3 mm, tandis que pour le dipôle radial des dipôles de faibles moments étaient retrouvés avec des valeurs élevées du critère à minimiser, et des erreurs de localisation variant entre 16.9 mm et 32.5 mm. Pour le modèle BEM où le cerveau était modélisé par une demi sphère, les dipôles tangentiels étaient localisés avec des erreurs de 6.5 mm et 9 mm, tandis que pour le dipôle radial l'erreur de localisation était de 17.9 mm avec un moment dipolaire d'amplitude égale à 20% de celle des dipôles tangentiels et une valeur du critère comparable à celles obtenues pour les dipôles tangentiels.

Tomita et al. [1996] ont placé 15 dipôles tangentiels dans un modèle réaliste à une couche et le champ magnétique engendré dans ce modèle a été calculé. Ces données ont ensuite été utilisées avec un modèle sphérique pour résoudre le problème inverse. Cette simulation a donné un biais de localisation croissant d'une moyenne de 2 mm pour les dipôles périphériques à 14 mm avec un maximum de 21 mm pour les dipôles les plus profonds. Le biais d'orientation était en moyenne de 5° pour les dipôles périphériques et de 18° avec un maximum de 63° pour les dipôles les plus profonds. Zanow [1997] a réalisé une simulation similaire. Il a montré que le biais de localisation variait de 1 mm pour des dipôles périphériques près du vertex à 25 mm pour des dipôles profonds dans le cortex frontal. Il a également résolu le problème inverse avec des modèles réalistes différents de celui utilisé pour le problème direct (obtenus en recommençant la segmentation des IRM), mais avec le même nombre de triangles. Les erreurs de localisation obtenues avec ces modèles réalistes mais de forme différente permettaient d'obtenir une amélioration de l'erreur de localisation par rapport au modèle sphérique dans presque tous les cas. Dans l'article de Meijs et al. [1988], c'est un modèle réaliste à quatre couches (scalp, os, liquide céphalo-rachidien (LCR), cerveau) qui est utilisé pour calculer le problème direct avec 20 dipôles tangentiels situés dans la région occipitale du cerveau. Le problème inverse est alors résolu soit avec un modèle sphérique, soit avec un modèle de quatre sphères excentriques, soit avec un modèle réaliste à une couche (le cerveau), soit avec un modèle réaliste à deux couches (cerveau et LCR). Les erreurs de
localisation et d'orientation moyennes obtenues sont de 4.5 mm et 3.2° pour le modèle sphérique, 3.8 mm et 3.2° pour le modèle sphérique excentrique, 4 mm et 3.7° pour le modèle réaliste à une couche et 3.7 mm et 3.3° pour le modèle réaliste à deux couches. Selon cette étude, en haut de la région occipitale il vaut mieux utiliser un modèle réaliste qu'un modèle sphérique ou à sphères excentrées, mais en bas de la région occipitale il vaut mieux utiliser un modèle sphérique.

## Évaluation des modèles réalistes sur un fantôme

Plusieurs études utilisent des fantômes physiques de la tête pour étudier la précision des modèles réalistes. Ces fantômes permettent de mieux tenir compte du bruit de mesure que les simulations, mais, tout comme les modèles simulés, ils n'ont pas nécessairement exactement les mêmes caractéristiques de conductivité que la tête d'un sujet vivant, et ne rendent pas compte du bruit physiologique.

Barth et al. [1986] ont placé 5 dipôles dans un fantôme constitué de la tête d'un cadavre humain dans laquelle le cerveau avait été remplacé par une gelée conductrice. Les erreurs de localisation obtenues avec un modèle sphérique variaient entre 1.3 mm pour les sources périphériques et 9.5 mm pour les sources profondes.

[Menninghaus et al., 1994] ont placé 4 dipôles dans un fantôme constitué d'un crâne en plastique non conducteur rempli de gelée conductrice. Les erreurs de localisation avec le modèle sphérique et variaient entre 3 mm et 9 mm pour 4 dipôles.

## Évaluation des modèles réalistes avec des dipôles implantés chez des patients épileptiques

Cohen et al. [1990] ont enregistré 12 dipôles implantés chez des patients épileptiques. Les erreurs de localisation variaient entre 4 mm et 40 mm. Toutefois 10 des dipôles avaient une orientation radiale et les données étaient très bruitées.

Balish et al. [1991] ont enregistré 19 dipôles implantés chez des patients épileptiques. Les erreurs de localisation variaient entre 3 mm et 40 mm. Les erreurs de localisation importantes correspondaient aux données les plus bruitées.

Rose et al. [1991] ont enregistré 3 dipôles implantés chez des patients épileptiques. Les erreurs de localisation variaient entre 1 mm et 4 mm.

## Carences des évaluations en MEG

Aucune étude aussi systématique que celle de Yvert et al. [1997] pour l'EEG ne semble avoir été réalisée en MEG. L'article de Meijs et al. [1988] montre que le nombre de couches dans un modèle réaliste pourrait avoir une importance qui n'existe pas dans le modèle sphérique. Contrairement aux modèles sphériques, dans les modèles réalistes, les composantes
radiales des dipôles ne sont pas muettes. Aucune étude ne semble avoir été faite sur l'effet de l'orientation des dipôles sur les erreurs de localisation et d'orientation avec les modèles réalistes.

1.3.2 Apport de l'interpolation linéaire du potentiel dans la méthode des intégrales du surface

Lorsque la méthode de Galerkin est utilisée (cf. section 1.2.3.2), on peut trouver une
formule de majoration de l'erreur commise par la méthode des intégrales de surface pour
notre problème [Dautray et Lions, 1985b] :

$$
\|V - \hat{U} \circ P^{-1}\|_{L^2(S)} \le C(h^{k+1} + h^m)\|V\|_{H^{k+1}(S)}
$$

Dans cette majoration, $V$ est le potentiel solution, $\hat{U}$ est la solution approchée obtenue par
la méthode des intégrales de surface, $P$ est la projection orthogonale qui relie localement la
surface réelle avec la surface approchée par des triangles, $h$ est la taille de la plus grande arête
du maillage, $k$ est le degré de la fonction d'interpolation polynômiale pour le potentiel
approché, et $m$ est le degré de la surface d'interpolation pour la surface approchée. Pour la
méthode de collocation, nous n'avons pu trouver aucune formule, mais il est vraisemblable
que les degrés d'interpolation du potentiel et de la géométrie $k$ et $m$ sont également des
facteurs cruciaux.

De nombreux auteurs [De Munck, 1992; Ferguson et al., 1994; Schlitt et al., 1995; Zanow, 1997] suggèrent qu'une interpolation linéaire ($k = 1, LINT$) pour le potentiel est préférable à l'approche centre de masse ($k = 0, COM$). En effet le nombre de sommets dans un maillage fermé est environ deux fois moins important que le nombre de triangles, ce qui permet une diminution du nombre d'inconnues, et donc du temps de calcul. De Munck [1992] propose une formule analytique pour le calcul de la matrice de la méthode des intégrales de surface dans le cas de l'interpolation linéaire (voir annexe C). Cependant cette formule ne s'applique qu'à des triangles plats ($m = 1$). Il propose d'étendre sa formule à des triangles courbes ($m = 2$) en décomposant chaque triangle courbe en plusieurs triangles plats, mais il s'agit d'une méthode approchée, et on perd alors les avantages de la formule analytique.

Une étude en EEG portant sur des patients épileptiques a montré qu'avec l'approche
interpolation linéaire ($m = 1, k = 1$) on obtient des erreurs de localisation de l'ordre de 8
millimètres en moyenne [Homma et al., 1994]. Cette étude a été réalisée en soumettant les
patients à un courant de faible intensité entre deux électrodes intracrâniennes pour réaliser un
dipôle artificiel. Elle présente l'avantage de tenir compte également de toutes les autres
sources d'erreur intervenant dans ce problème : erreur de placement des électrodes, erreur due
au modèle, au bruit.
Schlitt et al. [1995] ont mené une étude comparative par simulation entre l'approche centre de masse ($m = 1, k = 0$) et l'approche interpolation linéaire ($m = 1, k = 1$). Il résulte de cette étude réalisée avec huit dipôles tangentiels que l'approche centre de masse est plus précise que l'approche interpolation linéaire en ce qui concerne le potentiel mesuré sur le scalp. Cependant il semble que cette étude ait été réalisée sur des maillages obtenus par homothétie à partir du même maillage, or Yvert et al. [1995] indique que l'utilisation de trois maillages sphériques homothétiques conduit à des erreurs significativement moins importantes. Comme les maillages réalistes ne remplissent pas ces conditions de symétrie, il est nécessaire de s'affranchir de ce biais, par exemple en imposant des rotations différentes aux maillages de chaque couche.

Bien que cette interpolation semble, depuis la publication de De Munck [1992], avoir la préférence de nombreux auteurs [Ferguson et al., 1994; Schlitt et al., 1995; Zanow, 1997], nous n'avons pu trouver aucune étude systématique démontrant l'apport de l'interpolation linéaire du potentiel dans la méthode de collocation.

### 1.3.3 Apport de l'affinage local des maillages

La précision des modèles réalistes se dégrade lorsque les sources de courant sont proches de la surface, et des erreurs de localisation pouvant dépasser 1 centimètre peuvent alors être constatées en EEG [Yvert et al., 1996]. On peut espérer atténuer cet effet en affinant le maillage à proximité de la source, ainsi que le suggèrent Meijs et al. [1987b]. On pourrait penser que cette approche est irréaliste, puisqu'on ne connaît pas la position de la source a priori. Cependant, en fonction du protocole expérimental utilisé, une première localisation grossière peut être effectuée, et des maillages affinés adaptés aux régions explorées peuvent être utilisés. Pour les données de type EEG, une étude systématique réalisée par Yvert et al. [1996] sur des modèles réalistes (collocation, approche centre de masse) et sphériques simulés avec 24 dipôles montre que la méthode d'affinage est efficace, et permet de ramener l'erreur de localisation pour les dipôles dont la profondeur est inférieure à 3 millimètres. Le maillage optimal semble être obtenu lorsque le rayon de la zone d'affinage est égal à 3 centimètres. Zanow [1997] montre que l'affinage local des maillages apporte également une amélioration de la précision du problème direct dans le cas des données de type MEG, avec un modèle sphérique et un modèle réaliste à une couche (collocation, approche interpolation linéaire). Cependant, la stratégie d'affinage utilisée n'est pas du tout identique. Elle prend en considération la distribution du potentiel électrique pour réaliser le maillage (maillage adaptatif). Cette stratégie suppose de recalculer la matrice A de l'équation 1.2.1.12 et son inverse à chaque itération du problème inverse, ce qui demande beaucoup de temps de calcul.
C'est probablement la raison pour laquelle l'auteur ne présente pas de validation en problème inverse.

### 1.3.4 Apport du couplage MEG/EEG

Théoriquement, dans le modèle sphérique, les sources radiales et les composantes radiales des sources ne peuvent être enregistrées en MEG. De plus, l'utilisation d'un tel modèle conduit à des erreurs de localisation importantes lorsque les sources sont profondes pour la MEG [Barth et al., 1986; Roth et al., 1993; Menninghaus et al., 1994]. Une étude plus récente montre que si ces effets sont nettement moins importants lorsque l'on utilise des modèles réalistes, ils restent présents [Menninghaus et Lütkenhöner, 1995].

L'EEG enregistre toutes les sources, qu'elles soient radiales ou tangentielles, profondes ou périphériques, mais peut induire des erreurs de localisation plus importantes [Stok, 1987]. Stok montre cependant que ces deux types de données sont complémentaires, et que la résolution du problème inverse sphérique à 3 couches sur les deux modalités simultanément conduit à de bien meilleures localisations que lorsqu'une seule des deux modalités est utilisée. Cette minimisation simultanée est réalisée avec un critère à minimiser pondérant chaque voie par l'inverse de la variance du bruit sur cette voie. Les potentiels et champs magnétiques utilisés dans cette simulation ont été calculés analytiquement pour un modèle sphérique en ajoutant sur chaque voie un bruit gaussien.

Ce résultat est confirmé par l'étude de Mosher et al. [1993] pour un modèle sphérique à 4 couches. Ils montrent que la variance de l'erreur de localisation est plus faible et plus indépendante de l'orientation du dipôle lorsque les deux modalités sont utilisées simultanément. Cette étude porte sur un grand nombre de dipôles situés tous les millimètres sur un quart de plan dans une sphère de rayon 79 millimètres et avec des orientations dans le plan tangentiel espacées de 1° entre 0 et 179°. Cette fois encore, les potentiels et champs magnétiques de départ sont calculés analytiquement pour un modèle sphérique en ajoutant sur chaque voie un bruit gaussien, et les deux modalités sont ramenées à des grandeurs comparables en multipliant l'une d'elles par le rapport des variances des bruits.

Ces différents travaux montrent que le problème du couplage peut être abordé de différentes manières. Toujours dans le cas du problème inverse sphérique, Cohen et Cuffin [1987], et Hasson [1991] proposent trois méthodes différentes pour résoudre le problème inverse en utilisant simultanément les deux modalités. Cette fois encore, la validation de la méthode porte sur une simulation informatique avec des signaux bruités.

Dans sa thèse, Wieringa [1993] propose également une approche couplant ces deux modalités avec l'imagerie par résonance magnétique et la tomographie par émission de positons. Un exemple clinique est présenté. Cette fois, les différentes modalités sont utilisées
séparément pour le problème inverse sphérique. Cette étude portant sur un seul sujet montre que les différentes données apportent des informations complémentaires qu'il est utile de mettre en correspondance.

Une autre étude clinique a également été réalisée en utilisant les deux modalités [Buchner et al., 1994]. Dans cette étude, les minimisations ont été réalisées séparément, mais les données de type MEG ont permis de poser des contraintes pour la position et l'orientation des sources dans le cas de sources multiples. Avec ces contraintes, la solution au problème inverse pour l'EEG a pu être stabilisée. Cette fois encore, le modèle utilisé était un modèle sphérique. Le recalage des dipôles ainsi trouvés sur les images par résonance magnétique du cerveau a montré que ces résultats étaient en accord avec l'organisation somatotopique du sillon central.

Au début de cette thèse, aucune étude à notre connaissance n'utilisait les deux modalités simultanément, que ce soit avec le modèle sphérique ou avec un modèle réaliste, dans le cas de signaux réels ou obtenus par modélisation réaliste. Dans le cours de cette thèse, de nouveaux travaux ont été publiés, une bibliographie plus détaillée sera effectuée au chapitre 4. De nombreuses questions restent ouvertes dans ce domaine prometteur : comment coupler des signaux mesurés dans des unités physiques différentes (le microVolt et le femtoTesla), soumis à des bruits différents ? Quelle proportion d'électrodes doit-on utiliser par rapport au nombre de voies pour la MEG ? Quelle amélioration peut-on attendre du couplage en fonction de la position et de l'orientation des sources ? L'existence d'appareils permettant la mesure simultanée des deux modalités confirme la nécessité d'approfondir cette voie.

## 1.4 Plan de travail

Dans cette thèse, plusieurs hypothèses restrictives ont été choisies :

* Le modèle de sources est le modèle dipolaire. Actuellement, différents modèles de sources distribuées sont utilisés [Hämäläinen, 1992; Mosher et al., 1992; Wang et al., 1992; Ilmoniemi, 1993; Pascual-Marqui et al., 1994; Baillet et Garnero, 1997; Huang et al., 1997; Phillips et al., 1997]. Toutefois les modèles de sources dipolaires ponctuelles restent très utilisés. De plus, les méthodes que nous pouvons proposer pour améliorer la localisation avec un modèle dipolaire devraient s'appliquer également dans le cas des sources distribuées, alors que la réciproque n'est pas toujours vraie. Nous avons fait le choix du modèle dipolaire ponctuel afin que les résultats soient plus transposables. De plus, les techniques de sources distribuées sont encore en cours d'évolution.

* Les modèles réalistes sont des modèles à trois couches de conductivités homogènes, et la méthode numérique choisie pour résoudre le problème direct est la méthode des
intégrales de surface. Concernant les méthodes numériques, il y a principalement deux autres méthodes qui pourraient être utilisées. La méthode des différences finies est attirante par la facilité de la construction des maillages à partir des IRM. Toutefois, elle demande des temps de calcul beaucoup trop longs pour être utilisable en pratique. La méthode des éléments finis permet la prise en compte de l'anisotropie de l'os. Toutefois la génération de maillages volumiques des milieux conducteurs de la tête pose plus de problèmes que la génération de maillages surfaciques, et les valeurs des conductivités sont encore mal connues.

* Les aspects temporels ne sont pas étudiés. L'activation d'une source neuronale a toujours un décours temporel. L'exploitation des aspects temporels des signaux MEG et EEG permet donc une meilleure discrimination des différentes sources actives. Toutefois les améliorations apportées dans le cas statique s'ajoutent à celles apportées par le traitement temporel des signaux.

Ainsi que nous venons de le voir, les modèles réalistes semblent améliorer la localisation des sources de l'activité électrique cérébrale mesurée par EEG. Toutefois, un certain nombre de techniques pourraient permettre d'améliorer encore la précision des modèles réalistes. Ainsi, l'utilisation d'une interpolation plus précise du potentiel semble prometteuse. Au chapitre 2, l'apport de l'utilisation d'une interpolation linéaire du potentiel par rapport à une interpolation constante sera évalué.

Dans le cas de la MEG, l'influence de la géométrie du modèle et des différentes conductivités est réputée moins importante qu'en EEG, toutefois aucune évaluation systématique de l'erreur commise par le modèle sphérique lors de la résolution du problème inverse ne semble avoir été publiée. De plus, les modèles réalistes en MEG peuvent théoriquement détecter les composantes radiales des dipôles. C'est pourquoi l'apport des modèles réalistes et l'effet de l'orientation du dipôle en MEG sera étudié au chapitre 3.

Enfin le couplage des deux modalités MEG et EEG devrait permettre d'améliorer encore la précision de la reconstruction des sources de l'activité électrique cérébrale. Or toutes les études sur l'apport du couplage que nous avons pu trouver dans la littérature utilisent un modèle sphérique. Le chapitre 4 exposera les méthodes et les résultats obtenus dans des simulations par couplage MEG/EEG utilisant des modèles à géométrie réaliste.
<figure><img src="image_14.png" /></figure>
# Chapitre 2

## Amélioration de la méthode BEM en EEG
par une interpolation plus précise
## 2.1 L'interpolation linéaire du potentiel

Dans ce chapitre, nous comparons l'approche centre de masse avec l'approche interpolation linéaire pour la localisation des sources de courant intracérébrales à partir de la mesure du potentiel sur le scalp. Nous avons utilisé 56 dipôles sur les axes pariétal et temporal dans un modèle sphérique à trois couches, et nous avons comparé les solutions obtenues par la méthode des intégrales de surface avec la solution analytique.

Nous avons utilisé des maillages de sphères ayant subi des rotations aléatoires entre surfaces pour s'affranchir du cas particulier des maillages sphériques homothétiques donnant lieu à des erreurs artificiellement plus faibles [Yvert et al., 1995].

L'influence du nombre d'électrodes a également été étudiée. Nous avons trouvé que l'interpolation linéaire n'apporte pas de précision supplémentaire lorsqu'on se situe à modèle identique. Cependant, pour une même taille de matrice, l'interpolation linéaire permet une petite amélioration de l'erreur de localisation.

### 2.1.1 Introduction

Pour résoudre le problème direct dans le cas de modèles réalistes avec la méthode des intégrales de surface, la première approche qui ait été utilisée dans la littérature pour la méthode des intégrales de surface est l'approche centre de masse (COM, pour Center Of Mass). C'est l'approche la plus simple. Dans cette approche, le potentiel est supposé constant sur chaque triangle du maillage. Les points de discrétisation sont alors les centres de gravité des triangles. Cependant, dans un maillage d'une surface fermée, il y a en moyenne deux fois plus de triangles que de sommets de triangles. C'est pourquoi une approche interpolation linéaire (LINT, pour Linear

INTERpolation) semble plus avantageuse : si on suppose que le potentiel varie linéairement sur chaque triangle, il suffit de connaître ses valeurs aux sommets des triangles. Ainsi, pour une meilleure interpolation du potentiel, on obtient un système linéaire plus petit, donc plus rapide à résoudre, et aboutissant a priori à une meilleure précision.

Dans un article, De Munck montre que le calcul de la matrice du système linéaire peut se faire à l'aide d'une formule analytique [De Munck, 1992]. Cependant, à notre connaissance, seuls Schlitt et al. [1995] et plus récemment (postérieurement aux résultats présentés ici) Mosher et al. [1999] ont testé l'efficacité de l'approche LINT et l'ont comparée à l'approche COM. Schlitt et al. [1995] montrent que pour un petit nombre de dipôles, l'approche LINT donne de meilleurs résultats que l'approche COM pour estimer le potentiel sur la surface interne de l'os. En revanche, l'approche LINT donne de moins bons résultats que l'approche COM sur la surface du scalp. Ce test a cependant été réalisé sur des maillages de sphères
obtenus homothétiquement à partir du même maillage. Or, comme le soulignent Yvert et al. [1995], les maillages homothétiques conduisent à une erreur artificiellement plus faible pour la méthode COM. Mosher et al. [1999] ont également montré que l'approche LINT donne de moins bons résultats que l'approche COM pour estimer le potentiel sur le scalp, mais ils n'indiquent pas s'ils ont utilisé des maillages homothétiques ou non. Ferguson et Stroink [1997] trouvent des résultats comparables pour les erreurs obtenues avec les approches LINT et COM avec des maillages à une couche sphériques, ovoïdes ou cubiques. Aucun de ces articles ne propose d'évaluation de l'approche LINT en comparant les erreurs de localisation obtenues par résolution du problème inverse.

Nous allons présenter ici les résultats que nous avons trouvé lorsque nous avons testé l'approche interpolation linéaire en problème direct et en problème inverse. Bien que l'analyse des erreurs soit rendue difficile par la non-linéarité de sa dépendance aux différents paramètres (position du dipôle, orientation, nombre et position des électrodes), nous verrons que l'approche interpolation linéaire n'apporte que peu d'amélioration à la précision de localisation. Nous proposerons alors quelques pistes pour expliquer ce phénomène.

### 2.1.2 Méthode

Les équations matricielles que nous avons utilisées sont celles qui sont présentées dans l'article de De Munck [1992] (voir la section 1.2.3.2 et l'annexe C). Nous avons utilisé la technique de déflation de Lynn et Timlake [1968] pour rendre le système inversible et compenser le fait que le potentiel n'est défini qu'à une constante additive près. Nous avons également utilisé l'approche problème isolé suggérée par Hämäläinen et Sarvas [1989] et Meijs et al. [1989] après avoir vérifié qu'elle s'appliquait toujours dans le cas de l'interpolation linéaire du potentiel. Les développements mathématiques de cette approche sont décrits plus en détail en annexe de ce rapport.

Les évaluations ont été réalisées sur des maillages de surfaces sphériques de manière à pouvoir comparer le potentiel obtenu pour la méthode des intégrales de surface à celui donné par la solution analytique. Les conductivités sont supposées constantes à l'intérieur de chaque couche, et ont pour valeurs 0.45 S/m pour le scalp et le cerveau, et 0.005625 S/m (80 fois plus faible) pour l'os. Les couches sont délimitées par trois sphères concentriques de rayons 90 mm pour le scalp, 82.8 mm pour l'os et 78.3 mm pour le cerveau.

La méthode utilisée pour construire et caractériser les maillages est la même que celle utilisée dans l'article de Yvert et al. [1995]. Ils sont construits à l'aide de l'algorithme itératif de Delaunay, qui permet d'obtenir un maillage d'une surface monovaluée initialement représentée par un ensemble de points [Bertrand et al., 1991]. Ce maillage peut être affiné localement. Deux paramètres sont utilisés pour caractériser les maillages : la densité globale
(GD) et la densité locale (LD) de triangles par centimètres carrés respectivement sur toute la surface ou localement sur la surface affinée. La surface est affinée sur un angle d'ouverture de 40° (voir figure 2.1). Le nombre total de triangles par surface sera noté N. Les méthodes LINT et COM ont été testées à la fois pour le problème direct et pour le problème inverse. Le tableau 2.1 présente les maillages utilisés pour le problème direct, et le tableau 2.2 présente ceux utilisés pour le problème inverse.

<figure><img src="image_15.png" /><figcaption>FIG. 2.1: Exemple de maillage de sphère localement affiné sur un angle d'ouverture de 40°</figcaption></figure>

### 2.1.2.1 Évaluation du problème direct

L'erreur commise par la méthode des intégrales de frontière dans le problème direct dépend de plusieurs paramètres : la position et l'orientation du dipôle, la densité globale et la densité locale du maillage. Nous avons donc réalisé nos simulations pour différentes valeurs de ces paramètres. L'erreur commise sera caractérisée par deux critères introduits par Meijs et al. [1989] :

$$MAG = \sqrt{\frac{\int_{scalp} V_{num}^2}{\int_{scalp} V_{anal}^2}}$$

$$RDM = \sqrt{\int_{scalp} \left( \frac{V_{anal}}{\int_{scalp} V_{anal}^2} - \frac{V_{num}}{\int_{scalp} V_{num}^2} \right)^2}$$

Le MAG (pour MAGnification error) représente l'erreur d'amplitude, et le RDM (pour Relative Difference Measure) représente l'erreur topographique. Ces deux critères ont été calculés pour des dipôles dont la distance au centre varie entre 0 et 78 mm et pour des orientations du moment dipolaire radiales ou tangentielles par rapport à la surface du scalp. Ces dipôles sont situés sur l'axe z, et c'est autour de l'axe z que les maillages ont été
<table><caption>TAB. 2.1: Maillages utilisés pour le problème direct. GD = densité globale (triangles/cm²). LD = densité locale (triangles/cm²). N = nombre de triangles.</caption><thead><tr><th>Nom</th><th>GD</th><th>LD</th><th>N</th></tr></thead><tbody><tr><td>156elts</td><td>0.15</td><td></td><td>156</td></tr><tr><td>248elts</td><td>0.25</td><td></td><td>248</td></tr><tr><td>354elts</td><td>0.35</td><td></td><td>354</td></tr><tr><td>556elts</td><td>0.56</td><td></td><td>556</td></tr><tr><td>810elts</td><td>0.81</td><td></td><td>810</td></tr><tr><td>1328elts</td><td>1.33</td><td></td><td>1328</td></tr><tr><td>2334elts</td><td>2.33</td><td></td><td>2334</td></tr><tr><td>raf638elts</td><td>0.51</td><td>2.8</td><td>638</td></tr><tr><td>raf774elts</td><td>0.51</td><td>6.3</td><td>774</td></tr><tr><td>raf892elts</td><td>0.51</td><td>9.9</td><td>892</td></tr><tr><td>raf1112</td><td>0.51</td><td>17.9</td><td>1112</td></tr></tbody></table>

<table><caption>TAB. 2.2: Maillages utilisés pour le problème inverse. GD = densité globale (triangles/cm²). LD = densité locale (triangles/cm²). N = nombre de triangles.</caption><thead><tr><th>Nom</th><th>GD</th><th>LD</th><th>N</th></tr></thead><tbody><tr><td>604elts</td><td>0.6</td><td></td><td>604</td></tr><tr><td>604-900elts</td><td>0.6</td><td>6</td><td>900</td></tr><tr><td>900elts</td><td>0.9</td><td></td><td>900</td></tr><tr><td>900-1200elts</td><td>0.9</td><td>6.1</td><td>1200</td></tr><tr><td>1200elts</td><td>1.2</td><td></td><td>1200</td></tr></tbody></table>
localement affinés. Pour chaque dipôle et pour chaque maillage, les deux approches : COM et LINT ont été testées. Dans [Yvert et al., 1995], il a été montré que dans l'approche COM, l'utilisation de maillages sphériques dont les trois couches sont obtenues homothétiquement à partir du même maillage provoquait une diminution significative de l'erreur. En pratique, lors de l'utilisation de modèles réalistes, les triangles d'une surface sont orientés de manière aléatoire par rapport à ceux des surfaces voisines. C'est pourquoi nous donnerons des moyennes du MAG et du RDM pour 12 réalisations de trois rotations aléatoires différentes pour les trois sphères de chaque modèle.

## 2.1.2.2 Évaluation du problème inverse

Pour l'évaluation en problème inverse, trois montages d'électrodes selon le système international 10-20 avec 19, 32 et 63 électrodes ont été envisagés dans chaque cas. Ces montages sont représentés en figure 2.2. Pour chaque dipôle, le potentiel aux électrodes a été calculé avec la solution analytique. L'algorithme utilisé pour résoudre le problème inverse est le même que celui qui est utilisé dans [Yvert et al., 1995]. C'est une combinaison de l'algorithme de Marquardt avec une minimisation linéaire (voir section 1.2.4.2). La solution obtenue avec cet algorithme dépend fortement des valeurs initiales données à chaque paramètre. Pour nous affranchir de ce problème, et afin d'estimer uniquement les erreurs liées à l'approximation numérique ou au nombre d'électrodes, nous avons initialisé les paramètres du dipôle avec les paramètres réels du dipôle considéré.

<figure><img src="image_16.png" /><figcaption>FIG. 2.2: Position des électrodes sur le scalp.</figcaption></figure>

La simulation permet de calculer l'erreur de localisation qui est définie comme la distance entre les positions du dipôle considéré et du dipôle obtenu par l'algorithme de résolution du problème inverse. Cette erreur est calculée pour des dipôles dont la distance au centre varie entre 0 et 78 mm, et pour des positions sur l'axe pariétal et l'axe temporal représentés figure 2.3. Les distributions de potentiel sur l'axe pariétal sont mieux échantillonnées par les
configurations d'électrodes. Pour chaque axe, une rotation a été appliquée aux maillages localement affinés de manière à ce que l'axe se situe au centre de la zone affinée. Trois réalisations de trois rotations aléatoires différentes pour chaque couche ont été appliquées aux maillages afin d'éviter le cas particulier des maillages homothétiques.

<figure><img src="image_17.png" /><figcaption>FIG. 2.3: Positions des axes temporal et pariétal par rapport aux électrodes, en vue de derrière. Les cercles gras représentent le montage à 19 électrodes, tandis que les cercles fins représentent le montage à 32 électrodes.</figcaption></figure>

### 2.1.3 Résultats

#### 2.1.3.1 Problème direct

La figure 2.4 représente le RDM obtenu en moyenne pour 11 réalisations de rotations aléatoires différentes pour les trois sphères du modèle par rapport à ceux obtenus pour les modèles à maillages homothétiques. Les 7 modèles non affinés de la table 2.1 et 15 dipôles sont envisagés. On y voit que le cas homothétique dans l'approche COM donne de bien meilleurs résultats que les maillages ayant subi des rotations aléatoires, notamment pour les sources périphériques. Ce cas de symétrie ne peut être étendu aux maillages réalistes puisque le cerveau, l'os et le scalp n'ont pas des frontières homothétiques. Il convient donc de s'affranchir de ce biais pour la comparaison entre les approches COM et LINT en utilisant des maillages ayant subi des rotations aléatoires. Les résultats présentés ci-après correspondent donc à des moyennes de valeurs obtenues pour différentes réalisations de rotations aléatoires des maillages.

La figure 2.5 représente le MAG et le RDM en fonction de la distance des dipôles au centre du modèle, pour tous les dipôles et tous les maillages du tableau 2.1. Douze modèles et
<figure><img src="image_18.png" /><figcaption>FIG. 2.4: Comparaison entre les maillages homothétiques et les maillages ayant subi des rotations aléatoires des différentes couches pour l'approche COM.</figcaption></figure>

14 dipôles sont envisagés. Ainsi qu'on pouvait l'espérer, l'approche LINT donne un meilleur RDM que l'approche COM. Cependant, l'approche LINT donne également un MAG nettement moins bon que l'approche COM. Ce phénomène est mis en évidence sur la figure 2.6 qui représente le MAG et le RDM en moyenne pour tous les dipôles. On constate que tous les paramètres ont une influence sur l'erreur commise. L'influence de la taille du maillage sur le MAG est particulièrement plus prononcée dans le cas de l'approche LINT pour des maillages non affinés localement. Dans le cas des maillages affinés localement, les courbes sont plus regroupées pour les dipôles profonds, la densité globale des différents maillages étant identique. L'effet de l'affinage se fait sentir pour les dipôles distants de plus de 70 mm du centre.

Les comparaisons des figures 2.5 et 2.6 ont été réalisées en rapprochant les MAG et les RDM obtenus pour des modèles géométriques (maillages) identiques. Or les temps de calculs nécessaires à la résolution du problème direct dépendent de la taille de la matrice pour le système linéaire, qui est deux fois plus petite pour l'approche LINT. Une autre comparaison est donc effectuée pour une même taille de matrice à la figure 2.7.

Cette fois encore, le RDM est meilleur pour l'interpolation linéaire, en particulier pour les matrices de taille inférieure à 400, mais le MAG reste légèrement moins bon. Cependant, cette différence pour le MAG est très petite, et décroît rapidement lorsque la taille de la matrice augmente.
<figure><img src="image_19.png" /><figcaption>FIG. 2.5: Comparaison entre LINT et COM pour le problème direct en fonction de la distance au centre (en abscisse), de l'orientation du dipôle (radial ou tangentiel) et du modèle, localement affiné (L.A.) ou non (U.).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIG. 2.6: MAG et RDM en fonction du modèle pour différents dipôles : au centre, à 60 mm du centre et d'orientation radiale (rad 60 mm) ou tangentielle (tan 60 mm). Les carrés représentent l'erreur pour l'approche LINT, les ronds pour COM.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIG. 2.7: MAG et RDM en fonction de la taille de la matrice pour différents dipôles : au centre, à 60 mm du centre et d'orientation radiale (rad 60 mm) ou tangentielle (tan 60 mm). Les carrés représentent l'erreur pour l'approche LINT, les ronds pour COM.</figcaption></figure>

Cette comparaison effectuée sur le problème direct ne nous permet pas de conclure en faveur de l'approche LINT ou COM puisque les deux critères d'erreur donnent des résultats opposés. Il est donc nécessaire de mener une étude sur le problème inverse.
### 2.1.3.2 Problème inverse

Les figures 2.8 et 2.9 montrent l'évolution de l'erreur de localisation en moyenne pour les différentes réalisations de rotations aléatoires des maillages. Ces comparaisons, menées à modèle identique, montrent que l'interpolation linéaire détériore la précision de la localisation.

<figure><img src="image_22.png" /><figcaption>FIG. 2.8: Erreur de localisation en fonction du maillage utilisé, en moyenne pour 8 dipôles sur les axes temporal et pariétal, et pour trois rotations aléatoires des maillages. Trois montages de 19, 32 et 63 électrodes ont été utilisés. Les maillages 604-900 et 900-1200 sont des maillages affinés localement.</figcaption></figure>

Nous avons effectué un test statistique apparié non paramétrique de Wilcoxon pour comparer les méthodes COM et LINT pour toutes les réalisations de maillage et tous les dipôles. La méthode LINT donne, à maillages identiques, une erreur moyenne de 3.6 mm significativement plus grande que pour l'approche COM qui donne 3 mm (p < 0.0001).

La figure 2.9 montre que l'utilisation d'un plus grand nombre d'électrodes permet d'améliorer la localisation quelle que soit la méthode choisie. La figure 2.9 permet également une comparaison à taille de matrice identique. On voit alors que l'approche LINT donne des résultats légèrement meilleurs que l'approche COM. Ce résultat est confirmé par la figure 2.10, qui montre l'évolution de l'erreur de localisation en fonction de la taille de la matrice.

Le test statistique apparié non paramétrique de Wilcoxon pour les matrices de tailles 604 et 1200 montre que l'approche LINT donne cette fois de meilleurs résultats que l'approche COM avec des moyennes respectives de 3 mm et 3.3 mm (p = 0.0008).
<figure><img src="image_23.png" /><figcaption>FIG. 2.9: Erreur de localisation en fonction du nombre d'électrodes, en moyenne pour 8 dipôles sur les axes temporal et pariétal, et pour trois rotations aléatoires des maillages. Les modèles LINT 1200elts et LINT 2334elts ont respectivement à peu près la même taille de matrice que COM 604elts et COM 1200elts. Les modèles LINT 1200elts et COM 1200elts utilisent le même maillage.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIG. 2.10: Erreur de localisation en fonction de la taille de la matrice, en moyenne pour 8 dipôles sur les axes temporal et pariétal, et pour trois rotations aléatoires des maillages. Trois montages de 19, 32 et 63 électrodes ont été utilisés.</figcaption></figure>
### 2.1.4 Discussion

Cette simulation numérique aboutit à plusieurs résultats. Le premier d'entre eux est la confirmation de l'importance de ne pas garder pour l'évaluation des méthodes par simulation un modèle sphérique à trois couches dont les maillages sont obtenus homothétiquement à partir du même maillage. Cette observation nous amène à remettre en cause les résultats obtenus par Schlitt et al. [1995]. Les "RDM" (dont la définition n'est pas la même que la notre<fnref n="1" />) obtenus dans [Mosher et al., 1999] pour un temps de calcul identique laissent supposer qu'ils ont également utilisé des maillages homothétiques.

Par ailleurs, nous avons trouvé que l'approche LINT, pour un même modèle, donne un RDM un peu meilleur, mais un MAG moins bon, et aboutit à une erreur de localisation plus grande que l'approche COM. Cependant, si on effectue la comparaison *pour une même taille de matrice, ce qui correspond à un même temps de calcul, l'approche LINT améliore cette fois nettement le RDM, pour un MAG à peu près identique, et aboutit à une erreur de localisation légèrement plus faible*.

Comme le facteur limitant de la résolution du problème inverse est essentiellement le temps de calcul, on préférera donc l'approche LINT à l'approche COM.

Plusieurs explications peuvent être avancées pour expliquer le peu d'amélioration apportée par l'approche LINT, lorsqu'on se situe à modèle identique.

La première explication possible est que l'erreur peut dépendre non seulement du degré (k) du polynôme d'interpolation du potentiel, mais également du degré d'interpolation (m) de la surface. Dans l'approche LINT, nous n'avons pas augmenté le degré d'interpolation de la surface, puisque nous avons utilisé des triangles plats. Dans son article, De Munck [1992] fait allusion à ce problème. Il propose de le résoudre en découpant les triangles courbes en plusieurs triangles plats, de manière à pouvoir utiliser tout de même sa formule. Toutefois aucune évaluation de cette méthode n'est présentée. De plus cette méthode augmente le temps de calcul de la matrice.

Une autre explication qui peut être avancée pour expliquer ce phénomène est que l'erreur ne puisse être améliorée significativement en conservant une méthode de collocation. En effet, s'il existe [Dautray et Lions, 1985b] une formule de majoration de l'erreur commise par la méthode des intégrales de frontières avec des techniques variationnelles (Galerkin) par un $O(h^{k+1} + h^m)$ où $k$ est le degré des polynômes d'interpolation sur chaque élément de surface et $m$ est le degré d'interpolation géométrique des surfaces, nous n'avons pu trouver aucune formule semblable dans le cas de la méthode de collocation pour notre problème. L'utilisation de méthodes variationnelles a récemment été testée par Mosher et al. [1999], ils montrent une

<footnote n="1">¹Pour Mosher et al. [1999], $RDM = \sqrt{(b_{th} - b)^T (b_{th} - b) / (b_{th}^T b_{th})}$, où $b_{th}$ et $b$ représentent les valeurs trouvées analytiquement et numériquement à chaque voie pour l'EEG ou la MEG.</footnote>
amélioration de leur “RDM” par l’utilisation d’une interpolation linéaire par rapport à l’approche centre de masse. L’utilisation de méthodes variationnelles ne donne pas de meilleurs résultats que la collocation pour l’approche COM, mais donne de meilleurs résultats que la collocation pour l’approche LINT. Cependant l’apport des méthodes variationnelles a été évalué avec des maillages homothétiques, et en problème direct seulement.

## 2.2 Amélioration de l'interpolation géométrique

Nous avons testé ici l’approche proposée par De Munck [1992] en découpant des triangles courbes de degré 2 en plusieurs triangles plats. Les calculs sont expliqués dans l’annexe F. Les maillages de sphère avec les triangles plats et les triangles courbes sont représentés figure 2.11. Une comparaison entre l’approche LINT classique et cette approche est représentée figure 2.12.

<figure><img src="image_25.png" /><figcaption>FIG. 2.11: A gauche, un maillage de sphère avec 50 triangles plats. A droite, un maillage de sphère avec 50 triangles courbes. Chaque triangle courbe est subdivisé en 100 triangles plats. Certains triangles sont colorés de façon à pouvoir mieux les repérer pour les triangles courbes.</figcaption></figure>

La figure 2.12 montre que le fait de découper les triangles courbes en plusieurs triangles plats n'améliore pas sensiblement l'erreur de localisation, et l'aggrave même lorsque le maillage est grossier. Ceci est probablement dû au fait que l'on n'utilise plus véritablement une formule analytique pour calculer les coefficients de la matrice mais une formule approchée (voir annexe F).
<figure><img src="image_26.png" /><figcaption>FIG. 2.12: Erreur de localisation pour les triangles plats et pour les triangles courbes dans l'approche LINT. Deux maillages à 354 et 604 éléments par couche ont été utilisés, avec 63 électrodes. L'erreur est exprimée en fonction de la distance au centre pour les dipôles sur l'axe pariétal (pariet.) et l'axe temporal (tempor.), dans le cas d'un moment dipolaire orienté radialement (rad.) ou tangentiellement (tan.)</figcaption></figure>

### 2.3 Conclusions et perspectives

En conclusion, l'approche LINT, qui ne nécessite pas beaucoup de calculs supplémentaires et qui permet d'améliorer légèrement l'erreur de localisation à temps de calcul identique, pourra être utilisée pour remplacer l'approche COM. Cependant, on peut espérer améliorer davantage encore la précision de la méthode en utilisant soit une définition plus précise de la surface (avec de véritables triangles courbes) soit une méthode variationnelle.
<figure><img src="image_27.png" /></figure>
# Chapitre 3

## Apport des modèles réalistes en MEG
# 3.1 Introduction

## 3.1.1 Revue bibliographique

Le modèle le plus largement utilisé pour résoudre le problème inverse en MEG est le modèle sphérique. Il y a plusieurs raisons à cela :

* C'est un modèle très simple.

* Il existe une formule analytique au problème direct [Sarvas, 1987](voir annexe A).

* Les calculs avec la formule analytique sont très rapides.

* Des logiciels du commerce comme BESA® permettent depuis longtemps d'utiliser ce modèle.

* Ce modèle a la réputation de permettre une précision de localisation satisfaisante : dans le modèle sphérique, le champ magnétique ne dépend pas des valeurs des conductivités des milieux, dans la réalité, on suppose donc que le champ magnétique n'est pas beaucoup affecté par ces conductivités.

Des modèles plus élaborés composés de sphères non concentriques ont été introduits par Meijs et Peters [1987]. Cuffin [1991] a montré que la différence entre ces modèles et le modèle sphérique concentrique est minime en terme d'erreur de localisation (0.3 mm). Toutefois un seul dipôle situé sur un axe de symétrie a été étudié. Cuffin [1990] a également évalué la précision du modèle sphérique à partir de données simulées calculées avec la méthode des intégrales de surface dans des modèles à trois couches ellipsoïdaux ou dans un modèle sphérique où le cerveau était représenté par une demi-sphère. Pour trois dipôles (1 position et 3 orientations), le problème inverse était résolu avec un modèle sphérique analytique. Avec les modèles BEM ellipsoïdaux et les dipôles tangentiels l'erreur de localisation variait entre 2.6 mm et 12.3 mm et l'erreur d'orientation entre 0° et 3.2°, tandis que pour le dipôle radial des dipôles de faibles moments étaient retrouvés avec des valeurs élevées du critère à minimiser, et des erreurs de localisation variant entre 16.9 mm et 32.5 mm, et des erreurs d'orientation entre 58° et 76°. Pour le modèle BEM où le cerveau était modélisé par une demi sphère, les dipôles tangentiels étaient localisés avec des erreurs de localisation de 6.5 mm et 9 mm et des erreurs d'orientation de 0° et 6°, tandis que pour le dipôle radial les erreurs de localisation et d'orientation étaient de 17.9 mm et 73° avec un moment dipolaire d'amplitude égale à 20% de celle des dipôles tangentiels et une valeur du critère comparable à celles obtenues pour les dipôles tangentiels.

L'utilisation de modèles réalistes en MEG a été introduite par Hämäläinen [1987] et Meijs et al. [1987c]. Comme le champ magnétique est réputé peu dépendant des conductivités des milieux, on considère souvent qu'une seule couche (le cerveau) suffit pour modéliser la tête.
Ainsi dans [Zanow, 1997], une évaluation du biais de localisation introduit par le modèle sphérique par rapport à un modèle réaliste à une couche est réalisée. Il a montré que ce biais variait de 1 mm pour des dipôles périphériques près du vertex à 25 mm pour des dipôles profonds dans le cortex frontal. Dans [Tomita et al., 1996] une simulation similaire a donné un biais de localisation variant de 0.67 mm pour un dipôle occipital périphérique à 21 mm pour le dipôle frontal le plus profond. Le biais d'orientation de ces sources plutôt tangentielles variait entre 1° et 63°. Cependant, 12 des 15 dipôles étaient retrouvés avec une erreur d'orientation inférieure à 10°. Dans [Fuchs et al., 1998a], une même simulation a été effectuée, mais de façon beaucoup plus systématique (environ 3600 dipôles). Plusieurs configurations de capteurs MEG ont été étudiées : 31 gradiomètres axiaux situés soit sur la partie supérieure de la tête soit sur le côté, 109 gradiomètres axiaux (configuration tête entière), ou 122 gradiomètres planaires (tête entière). Pour chaque position de dipôle, les trois orientations propres ont été envisagées (voir section 1.2.2). Dans les orientations correspondant aux deux valeurs propres les plus grandes (orientations quasi-tangentielles), les erreurs de localisation variaient entre 0.2 mm et 70.9 mm (5.4 mm en moyenne) avec les 122 gradiomètres planaires, entre 0.2 mm et 77 mm (6.7 mm en moyenne) avec les 109 gradiomètres axiaux. Dans le cas des 31 gradiomètres axiaux, les erreurs de localisation les plus faibles étaient obtenues pour les sources situées à proximité des canaux, tandis que des erreurs de localisation très importantes étaient obtenues dans les régions non couvertes par les canaux. Dans l'orientation correspondant à la valeur propre la plus faible (orientation quasi-radiale) les erreurs de localisation étaient beaucoup plus importantes que dans les orientations quasi-tangentielles : en moyenne 44.1 mm pour les 122 gradiomètres planaires et 54 mm pour les 109 gradiomètres axiaux.

La précision du modèle sphérique n'a pas seulement été évaluée à partir de simulations. Des fantômes réalistes de la tête ont également été utilisés : dans [Barth et al., 1986], 5 dipôles ont été placés dans un fantôme constitué de la tête d'un cadavre humain dans laquelle le cerveau avait été remplacé par une gelée conductrice. Les erreurs de localisation obtenues avec un modèle sphérique et une voie MEG unique placée séquentiellement à 81 positions différentes variaient entre 1.3 mm pour les sources périphériques et 9.5 mm pour les sources profondes. Dans [Menninghaus et al., 1994] le fantôme était constitué d'un crâne en plastique non conducteur rempli de gelée conductrice. Il s'agissait donc d'un fantôme à une couche. Les erreurs de localisation avec le modèle sphérique et 37 gradiomètres axiaux variaient entre 3 mm et 9 mm pour 4 dipôles. Dans [Leahy et al., 1998], le fantôme était constitué d'un crâne humain imprégné et rempli de gelée conductrice recouvert d'une couche de latex afin d'obtenir un fantôme à trois couches. Les erreurs de localisation avec le modèle sphérique et l'algorithme R-MUSIC variaient entre 1 mm et 17 mm pour 32 dipôles. Les erreurs de localisation avec la méthode des intégrales de surface étaient améliorées de 0.5 mm
seulement. Toutefois le niveau de bruit dans leurs enregistrement était vraiment très important (leur rapport signal sur bruit en MEG était en moyenne de 2). De plus, leurs données étant disponibles sur Internet (http://neuroimage.usc.edu) nous avons pu constater qu'il y avait une intersection entre les couches des maillages, ce qui n'est pas autorisé dans la méthode des intégrales de surfaces.

Le modèle sphérique a également été évalué à partir de l'enregistrement de dipôles implantés créés entre deux dérivaions d'électrodes intracérébrales de profondeur chez des patients épileptiques [Cohen et al., 1990; Balish et al., 1991; Rose et al., 1991]. Des erreurs de localisation variant surtout en fonction de la profondeur du dipôle entre 3 mm et 20 mm, avec une moyenne selon les articles variant entre 8 mm et 10 mm ont été trouvées lorsque le rapport signal sur bruit était supérieur à 10.

Le tableau 3.1 résume les résultats trouvés dans la littérature sur l'erreur de localisation avec le modèle sphérique.

Dans le modèle sphérique, les composantes radiales des dipôles ne créent pas de champ magnétique et ne peuvent donc être estimées. Il est fréquemment admis que les modèles réalistes sont suffisamment proches du modèle sphérique pour que le champ magnétique créé par les composantes radiales des dipôles soit négligeable. C'est pourquoi peu d'études s'intéressent à l'erreur d'orientation des dipôles. Nous avons vu que Tomita et al. [1996] trouvaient des erreurs d'orientation inférieures à 10° pour 12 de ses 15 dipôles, toutefois tous ces dipôles ont une orientation plutôt tangentielle. Leahy et al. [1998] ont calculé uniquement la composante tangentielle de l'erreur d'orientation qui était en moyenne de 10° pour les 32 dipôles dans un fantôme à trois couches avec le modèle sphérique et de 8° avec un modèle réaliste. À partir de dipôles implantés dans le cerveau de lapins vivants, Melcher et Cohen [1988] ont trouvé une atténuation pour les dipôles radiaux de 0.17 ± 0.07 par rapport aux dipôles tangentiels. Menninghaus et Lütkenhöner [1995] ont également montré que les sources radiales pouvaient créer en réalité des champs magnétiques non négligeables, le rapport entre l'amplitude du champ créé par un dipôle tangentiel et par un dipôle radial variant de 8 pour les dipôles périphériques à 2 pour les dipôles profonds. De plus, ils ont montré un exemple de topographie de champ magnétique très différent pour un dipôle presque radial de ce que prédisait le modèle sphérique. Dans son modèle comportant un cerveau hémisphérique, Cuffin [1990] trouve que le champ magnétique engendré par un dipôle radial comporte des valeurs significatives (atténué de 22% par rapport à un dipôle tangentiel) et que la carte du champ magnétique ressemble à celle engendrée par un dipôle tangentiel.
<table><caption>TAB. 3.1: Erreurs de localisation reportées dans la littérature en utilisant le modèle sphérique. SNR : rapport signal sur bruit.</caption><thead><tr><th>Auteurs</th><th>Données MEG</th><th>Nombre de dipôles</th><th>Profondeur (mm) des dipôles</th><th>Orientation des dipôles</th><th>Érreur de localisation (mm)</th></tr></thead><tbody><tr><td>[Cuffin, 1991]</td><td>simulation, modèle de sphères non concentriques</td><td>1</td><td>31</td><td>Tangentiel</td><td>0.3</td></tr><tr><td>[Cuffin, 1990]</td><td>simulation, 5 modèles ellipsoïdaux à 3 couches</td><td>3</td><td>31</td><td>2 tangentiels, 1 radial</td><td>2.6-12.3, 16.9-32.5</td></tr><tr><td>[Cuffin, 1990]</td><td>simulation, scalp et os sphériques, cerveau hémisphérique</td><td>3</td><td>31</td><td>2 tangentiels, 1 radial</td><td>6.5, 9.0, 17.9</td></tr><tr><td>[Zanow, 1997]</td><td>simulation, modèle réaliste à une couche</td><td>30</td><td>7 - 53</td><td>Tangentiels</td><td>1 - 25</td></tr><tr><td>[Tomita et al., 1996]</td><td>simulation, modèle réaliste à une couche</td><td>15</td><td>≈ 15-55</td><td>Tangentiels</td><td>0.7 - 21.4</td></tr><tr><td>[Fuchs et al., 1998a]</td><td>simulation, modèle réaliste à une couche</td><td>≈ 3600</td><td>≈ 0-80</td><td>2/3 tangentiels, 1/3 radiaux</td><td>0.2-70.9, 13.1-87.9</td></tr><tr><td>[Barth et al., 1986]</td><td>Cadavre humain</td><td>5</td><td>≈ 20 - 75</td><td>Variable</td><td>1.3 - 9.5</td></tr><tr><td>[Menninghaus et al., 1994]</td><td>Fantôme à une couche</td><td>4</td><td>10 - 30</td><td>Tangentiels</td><td>3 - 9</td></tr><tr><td>[Leahy et al., 1998]</td><td>Fantôme à trois couches</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>Dipôles implantés</td><td>19</td><td>20 - 63</td><td>Tangentiels</td><td>3 - 40 3-20 si SNR &gt; 10</td></tr><tr><td>[Cohen et al., 1990]</td><td>Dipôles implantés</td><td>12</td><td>≈ 30 - 45</td><td>10 radiaux 2 tangentiels</td><td>4 - 40 4-13 si SNR &gt; 1</td></tr><tr><td>[Rose et al., 1991]</td><td>Dipôles implantés</td><td>3</td><td>≈ 15 - 45</td><td>Tangentiels</td><td>1 - 4</td></tr></tbody></table>
## 3.1.2 plan de travail

Dans [Yvert et al., 1997], une évaluation systématique du biais de localisation induit par l'utilisation d'un modèle sphérique en EEG a été réalisée. Nous présentons ici une étude similaire dans le cas de données MEG : le problème direct est résolu pour 2003 dipôles situés perpendiculairement à la surface corticale pour un casque MEG composé de 143 gradiomètres axiaux en utilisant un modèle réaliste à trois couches très finement maillé. Le problème inverse est ensuite résolu sur ces données simulées en utilisant soit un modèle sphérique, soit un modèle réaliste à une couche, soit un modèle réaliste à trois couches plus grossier que celui utilisé pour le problème direct. L'erreur de localisation, mais aussi l'erreur d'orientation et l'erreur sur l'amplitude du moment dipolaire sont présentées. Cette étude a fait l'objet d'un article [Crouzeix et al., 1999b](voir annexe G).

L'effet de l'orientation du dipôle sur les topographies de champ magnétique, sur l'amplitude du signal mesuré, et finalement sur les erreurs de localisation et d'orientation seront plus particulièrement approfondis dans un second temps.
## 3.2 Évaluation de l'apport des modèles réalistes pour la résolution du problème inverse

### 3.2.1 Méthode

#### 3.2.1.1 Génération des modèles réalistes

Les maillages des différentes surfaces : scalp, crâne et cerveau nécessaires à la méthode des intégrales de surface ont été créés en trois étapes :

1. Pour chaque image d'une pile d'IRM, les contours des différentes surfaces ont été extraits.

2. Les contours des différentes surfaces ont été reliés par un maillage très fin (environ 20000 nœuds) selon la procédure décrite dans [Yvert et al., 1995].

3. La technique de simplification de maillage proposée par Algorri et Schmitt [1996] a été utilisée pour réduire le nombre de triangles à une quantité utilisable en pratique.

Comme on utilise la méthode des intégrales de surface avec une interpolation linéaire du potentiel (LINT), le nombre de nœuds dans les maillages utilisés dans cette étude est un paramètre important, puisqu'il détermine la taille des équations matricielles à résoudre.

#### 3.2.1.2 Modèle utilisé pour créer les données simulées

Les données de champ magnétique simulées ont été obtenues en appliquant la méthode des intégrales de surface avec l'approche problème isolé [Meijs et al., 1987a; Hämäläinen, 1987] et une interpolation linéaire du potentiel (LINT) [De Munck, 1992; Ferguson et al., 1994]. Le maillage utilisé pour ces calculs (REF, voir table 3.2) était constitué d'environ 4500 nœuds pour une densité moyenne assez élevée de 3 triangles/cm². Les valeurs normalisées des conductivités pour les différentes couches étaient de 1,1/80 et 1 respectivement pour le scalp, le crâne et le cerveau, ce qui correspond aux valeurs suggérées dans [Rush et Driscoll, 1968; Homma et al., 1995].

#### 3.2.1.3 Configuration des capteurs MEG

Les données simulées ont été calculées avec une configuration de capteurs MEG formant un casque couvrant la tête entière et constituée de 143 gradiomètres axiaux de la société CTF. Ce casque a été placé de manière réaliste autour du modèle REF en conservant la distance entre le scalp et les capteurs aussi constante que possible. La distance minimum entre le scalp et les capteurs était de 19 mm. La figure 3.1 montre le positionnement des capteurs autour du modèle REF.
<figure><img src="image_28.png" /><figcaption>FIG. 3.1: Procédure de simulation. Les modèles utilisés pour le problème direct et le problème inverse sont représentés ainsi que les positions des spires collectrices.</figcaption></figure>
<table><caption>TAB. 3.2: Maillages utilisés pour le problème direct et le problème inverse</caption><thead><tr><th>Nom du modèle</th><th>Type</th><th>Nombre de couches</th><th>Nombre de nœuds<br>par couche</th><th>Densité de triangles<br>(tri./cm²)<br>scalp-crâne-cerveau</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center; font-weight: bold;">Problème direct</td></tr><tr><td>REF</td><td>Réaliste</td><td>3</td><td>1500</td><td>2.1 - 3 - 4</td></tr><tr><td colspan="5" style="text-align: center; font-weight: bold;">Problème inverse</td></tr><tr><td>GLSPH</td><td>Sphère globale approchant au mieux la surface du scalp</td><td>-</td><td>-</td><td>-</td></tr><tr><td>LOCSPH</td><td>Sphère locale approchant au mieux la surface du scalp</td><td>-</td><td>-</td><td>-</td></tr><tr><td>BRAINGLSPH</td><td>Sphère globale approchant au mieux la surface du cerveau</td><td>-</td><td>-</td><td>-</td></tr><tr><td>BRAINLOCSPH</td><td>Sphère locale approchant au mieux la surface du cerveau</td><td>-</td><td>-</td><td>-</td></tr><tr><td>1L900</td><td>Réaliste</td><td>1</td><td>900</td><td>2.6</td></tr><tr><td>3L300</td><td>Réaliste</td><td>3</td><td>300</td><td>0.5 - 0.7 - 1.2</td></tr><tr><td>3L900</td><td>Réaliste</td><td>3</td><td>900</td><td>1.9 - 1.8 - 2.6</td></tr></tbody></table>

### 3.2.1.4 Positions des dipôles

Les 2003 dipôles étaient utilisés dans [Yvert et al., 1997], ce qui permettra ultérieurement une comparaison entre les résultats obtenus en MEG avec les résultats obtenus en EEG. Leurs positions ont été tirées aléatoirement à l'intersection de la surface corticale et de 10 coupes horizontales séparées les unes des autres de 10 mm (voir figure 3.2). L'orientation des dipôles a été choisie perpendiculairement à la surface corticale pour chaque position. Les profondeurs des dipôles variaient entre 0.9 mm et 46.9 mm pour une moyenne de 16.4 mm et un écart-type de 9.8 mm.

L'amplitude du champ magnétique engendré par chaque dipôle est caractérisée par sa valeur RMS, représentée figure 3.3.

$$RMS = \sqrt{\frac{1}{143} \sum_{i=1}^{143} B_i^2}$$

Pour caractériser l'orientation des dipôles, nous avons défini un angle de tangentialité compris entre 0° et 90°. Il a été défini comme l'angle entre le moment dipolaire et sa composante radiale dans le modèle sphérique GLSPH défini ci-dessous. Par la suite, nous adopterons les dénominations suivantes :

* dipôle radial : dipôle dont l'angle de tangentialité est compris entre 0° et 5°.
<figure><img src="image_29.png" /><figcaption>FIG. 3.2: Une image 3D de la surface corticale segmentée a été superposée sur une coupe IRM sagittale. La projection des positions des dipôles est représentée par des points rouges.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIG. 3.3: Pour chaque position de dipôle, la valeur RMS correspondante est codée sous la forme d'un point coloré : rouge si $0 \text{ fT} < \text{RMS} < 0.5 \text{ fT}$; jaune si $0.5 \text{ fT} < \text{RMS} < 0.9 \text{ fT}$; vert si $0.9 \text{ fT} < \text{RMS} < 1.3 \text{ fT}$; bleu si $1.3 \text{ fT} < \text{RMS} < 14 \text{ fT}$. Dans cette figure, chaque couleur correspond à un quartile : il y a autant de points dans chaque couleur.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>FIG. 3.4: Pour chaque position de dipôle, l'angle de tangentialité correspondant est codé sous la forme d'un point coloré : rouge entre 0° et 30° (plutôt radial) ; jaune entre 30° et 60° ; vert entre 60 et 90° (plutôt tangentiel).</figcaption></figure>

* dipôle tangentiel : dipôle dont l'angle de tangentialité est compris entre 85° et 90°.

Parmi les 2003 dipôles, 63% avaient un angle de tangentialité compris entre 60° et 90° (voir figure 3.4).

### 3.2.1.5 Modèles utilisés pour le problème inverse

#### Modèles sphériques

Quatre modèles sphériques ont été utilisés pour résoudre le problème inverse (voir figure 3.1 et table 3.2) :

1. un modèle sphérique global approchant au mieux la surface du scalp couverte par le casque MEG (GLSPH).

2. un modèle sphérique semblable à GLSPH mais approchant au mieux la surface du cerveau couverte par le casque MEG (BRAINGLSPH).

3. un modèle sphérique local (LOCSPH) approchant au mieux la surface du scalp couverte par une sélection des voies MEG voisines de celles enregistrant le minimum et le maximum du signal. Cette surface a été déterminée comme suit :

* Soit $d$ la distance entre les deux spires mesurant le minimum et le maximum du signal.

* L'ensemble des spires dont la distance au milieu de ces deux spires était inférieure à $2 * d$ a été sélectionné. Si plus de la moitié des spires étaient sélectionnés, seul le modèle global était retenu, et le dipôle correspondant était écarté lors des comparaisons entre les différents modèles sphériques.
<figure><img src="image_32.png" /><figcaption>FIG. 3.5: Procédure adoptée pour choisir une sphère locale</figcaption></figure>

* La sphère locale approchait au mieux la surface du scalp dont la distance au milieu des 2 spires était inférieure à $2 * d$ (voir figure 3.5).

4. un modèle sphérique local semblable à LOCSPH mais approchant au mieux la surface du cerveau et non du scalp (BRAINLOCSPH).

Nous avons retenu cette procédure pour choisir une sphère locale car elle ne requiert pas la connaissance a priori de la position du dipôle et peut facilement être mise en œuvre en pratique lorsqu'une numérisation de la surface de la tête est disponible. Il est fréquent dans la pratique qu'aucune IRM ne soit disponible, et que seule une numérisation de la surface de la tête puisse être utilisée pour déterminer la sphère. Dans le cas contraire, on préfère généralement choisir une sphère approchant au mieux la surface du cerveau.

En plus de ces quatre modèles, l'influence de la position du centre de la sphère a été testée en utilisant six modèles sphériques translatis à partir de BRAINLOCSPH d'un centimètre dans les six directions axiales.

Dans le cas des modèles sphériques, l'orientation du dipôle est contrainte à rester tangentielle lors de la résolution du problème inverse.

## Modèle réaliste à une couche

Hämäläinen et Sarvas [1989] ont étudié les différences topographiques entre des solutions du problème direct obtenues pour un petit nombre de sources avec des modèles réalistes à une et trois couches. Ils sont partis de cette étude pour suggérer que, puisque la conductivité de l'os est très inférieure à celles du scalp et du cerveau, une seule couche suffisait pour un modèle réaliste de la tête en MEG. Nous avons utilisé un modèle à une couche (1L900) composé de 900 nœuds et ayant une densité de 2.5 triangles/cm² (voir figure 3.1 et table 3.2).
## Modèle réaliste à trois couches

Nous avons également utilisé deux modèles réalistes à trois couches : l'un plutôt grossier (3L300) et l'autre plus fin (3L900) avec respectivement 300 et 900 nœuds par couche (voir figure 3.1 et table 3.2). Les modèles 1L900 et 3L300 correspondent à des tailles de matrices et donc des temps de résolution identiques pour le problème direct.

### 3.2.1.6 Résolution du problème inverse

La résolution du problème inverse consiste à trouver les paramètres de position et d'orientation d du dipôle qui minimisent le critère suivant :

$$S(\mathbf{d}) = \sum_{i=1}^{N} (B_{REF}^{i} - B_{modelINV}^{i}(\mathbf{d}))^2$$

où :

* N est le nombre de voies MEG

* $B_{REF}^i$ est le champ magnétique calculé avec le modèle REF à la voie i.

* $B_{modelINV}^i(\mathbf{d})$ est le champ magnétique produit à la voie i par un dipôle de paramètres d pour un modèle inverse donné.

Pour cette minimisation, nous avons utilisé l'algorithme mixte comportant une partie linéaire et une partie non linéaire décrit à la section 1.2.4.2. Les solutions obtenues par cet algorithme dépendent fortement des valeurs initiales données aux paramètres des dipôles. Pour nous affranchir de ce problème, nous avons initialisé ces paramètres avec les paramètres des dipôles originaux. Pour 186 dipôles périphériques, le maillage 3L300 était tellement grossier que ces dipôles tombaient en dehors du maillage de la surface du cerveau. Ils ont donc été exclus dans les analyses statistiques. Par ailleurs, pour 50 dipôles, l'erreur résiduelle $S(\mathbf{d})$ divisée par le carré de la valeur RMS du signal restait supérieure à 30% à l'issue du problème inverse, ce qui signifie qu'aucune solution raisonnable n'avait été trouvée. Ces cas, qui arrivaient surtout avec les modèles sphériques, ont également été exclus des analyses statistiques.

### 3.2.1.7 Définitions des critères d'erreur

* L'erreur de localisation ($D$) est la distance en millimètres entre le dipôle reconstruit et le dipôle d'origine.

* L'erreur d'orientation ($\theta$) est l'angle en degrés entre les deux moments dipolaires.

* L'erreur sur l'amplitude du moment ($\Delta M$) est la différence entre l'amplitude du moment dipolaire du dipôle reconstruit et celle du dipôle original, en pourcentage de
l'amplitude du dipôle original. Les valeurs négatives et positives de cette erreur correspondent respectivement à la sous-estimation et la surestimation de l'amplitude.

## 3.2.2 Résultats

### 3.2.2.1 Comparaison des différents modèles sphériques

Pour 824 dipôles (dont la profondeur était généralement inférieure à 2 cm) la sphère locale (LOCSPH) et la sphère globale (GLSPH) étaient différentes. Pour tous ces dipôles, les sphères approchant au mieux la surface cérébrale (BRAINGLSPH) et (BRAINLOCSPH) étaient également différentes. La distance entre les centres des sphères approchant au mieux le scalp variait entre 2.5 mm et 58 mm avec une moyenne de 21 mm et un écart-type de 11.5 mm. La distance entre les centres des sphères approchant au mieux le cerveau variait entre 2.7 mm et 39 mm avec une moyenne de 14.7 mm et un écart-type de 6.9 mm. Pour la comparaison entre GLSPH, BRAINGLSPH, LOCSPH et BRAINLOCSPH, seuls ces 824 dipôles ont été pris en compte. Les erreurs obtenues avec ces modèles sont représentées figure 3.6.

#### Erreurs de localisation

dans la partie supérieure de la tête (coupes 1 à 5 de la figure 3.2), les erreurs de localisation obtenues en moyenne pour GLSPH, LOCSPH, BRAINGLSPH et BRAINLOCSPH étaient respectivement de 3.5 mm, 2.4 mm, 3.1 mm et 2.5 mm. La différence de 1.1 mm entre GLSPH et LOCSPH était significative (test t apparié : *p* < 0.0001), de même que la différence de 0.6 mm entre BRAINLOCSPH et BRAINGLSPH. En revanche, la différence de 0.1 mm entre LOCSPH et BRAINLOCSPH ne l'était pas (*p* = 0.29).

Dans la partie inférieure de la tête (coupes 6 à 10), les erreurs de localisation étaient plus importantes : 5.3 mm, 8.7 mm, 4.6 mm et 4 mm, respectivement pour GLSPH, LOCSPH, BRAINGLSPH, et BRAINLOCSPH. Cet effet vertical était significatif (test t non apparié : *p* < 0.0001). Une analyse plus détaillée dans cette zone montre que :

* Dans la partie arrière du cerveau, tous les modèles sphériques à l'exception de BRAINLOCSPH donnent une erreur de localisation moyenne entre 3.7 mm et 3.9 mm. Il n'y a pas de différence significative entre ces modèles (*p* > 0.6). BRAINLOCSPH donne une erreur de localisation moyenne de 2.9 mm ce qui est statistiquement meilleur que les autres modèles (*p* < 0.0001).

* Dans la partie antérieure du cerveau, les erreurs de localisation pour GLSPH, LOCSPH, BRAINGLSPH, et BRAINLOCSPH étaient respectivement de 6.5 mm, 13.4 mm, 5.2 mm et 5 mm. Toutes ces erreurs étaient significativement différentes (*p* < 0.0001) sauf entre BRAINLOCSPH et BRAINGLSPH (*p* = 0.62).
A

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

B

<figure><img src="image_34.png" /><figcaption>FIG. 3.6: Erreurs de localisation, d'orientation et d'amplitude du moment dipolaire en moyenne en fonction du numéro de coupe (A) et de la tangentialité (B). Les barres d'erreurs correspondent à une erreur standard. Seuls les 824 dipôles pour lesquels les sphères locales et globales diffèrent et l'erreur résiduelle est inférieure à 30% sont représentés.</figcaption></figure>
Nous avons fait des translations d'un centimètre de BRAINLOCSPH dans 6 directions parallèles aux axes du repère. L'erreur de localisation était affectée en moyenne de 0.8 mm. Aucune des 6 directions ne donnait de résultats significativement meilleurs.

## Erreurs d'orientation

Les différents modèles sphériques donnaient à peu près les mêmes erreurs d'orientation (figure 3.6). La position verticale du dipôle ne semblait pas avoir d'effet particulier sur cette erreur. Notons que chaque coupe contenait des dipôles d'orientations différentes. Ainsi que l'on pouvait s'y attendre, l'erreur d'orientation pour les dipôles radiaux était de 90° puisque le dipôle était contraint à rester tangentiel pendant la minimisation. Pour les dipôles tangentiels, l'erreur d'orientation descendait à 10°.

## Erreur sur l'amplitude du moment dipolaire

Ainsi que le montre la figure 3.6, l'erreur sur l'amplitude du moment dipolaire était affectée par l'erreur sur l'orientation du dipôle. Ainsi pour les dipôles radiaux, l'erreur sur l'amplitude du moment dipolaire était de -100%, et elle tombait entre -15% et 15% pour les dipôles tangentiels.

### 3.2.2.2 Erreurs obtenues avec les modèles réalistes

Les données statistiques indiquées dans cette section concernent les 1768 dipôles qui tombaient à l'intérieur de tous les maillages de cerveau et donnaient des erreurs résiduelles inférieures à 30%. Les résultats ont également été comparés à ceux obtenus avec BRAINLOCSPH et GLSPH.

## Erreurs de localisation

Les erreurs de localisation pour les modèles réalistes sont représentées figure 3.7, figure 3.8 et figure 3.10.

* Pour 1L900, l'erreur de localisation moyenne croissait de 2.4 mm dans la partie supérieure de la tête (coupes 1 à 5) à 4.9 mm dans la partie inférieure de la tête (coupes 6 à 10). Elle passait de 8.3 mm à 2.5 mm respectivement pour les dipôles radiaux et tangentiels.

* Pour 3L300 l'erreur de localisation moyenne croissait de 1.3 mm dans les parties supérieures et médianes de la tête (coupes 1 à 8) à 3.5 mm dans la partie inférieure de la tête (coupes 9 et 10). Elles variaient entre 4 mm et 0.97 mm respectivement pour les dipôles radiaux et tangentiels.
<figure><img src="image_35.png" /><figcaption>FIG. 3.7: Erreurs de localisation (A), d'orientation (B) et d'amplitude du moment dipolaire (C) en fonction de la tangentialité. Les barres d'erreurs correspondent à une erreur standard. Seuls les 1768 dipôles qui tombaient à l'intérieur de tous les maillages de cerveau et donnaient des erreurs résiduelles inférieures à 30% ont été pris en compte.</figcaption></figure>
* Pour 3L900 l'erreur de localisation moyenne croissait de 0.4 mm dans les parties supérieures et médianes de la tête (coupes 1 à 8) à 1.4 mm dans la partie inférieure de la tête (coupes 9 et 10). Elles variaient entre 2.6 mm et 0.4 mm respectivement pour les dipôles radiaux et tangentiels.

Toutes ces différences entre les parties supérieures et inférieures de la tête étaient statistiquement significatives ($p < 0.0001$). Les dipôles étaient localisés avec une erreur inférieure à 3 mm pour 59%, 90%, et 98% des cas avec respectivement 1L900, 3L300, et 3L900.

## Erreurs d'orientation

Les erreurs de localisation pour les modèles réalistes sont représentées figure 3.7, figure 3.9 et figure 3.10.

* Pour 1L900, l'erreur d'orientation moyenne variait avec la tangentialité entre 13° et 31°.

* Pour 3L300, l'erreur d'orientation moyenne variait avec la tangentialité entre 3° et 21°.

* Pour 3L900, l'erreur d'orientation moyenne variait avec la tangentialité entre 2° et 8°.

Les orientations des dipôles étaient retrouvées avec une erreur inférieure à 10° pour 27%, 60%, et 87% des dipôles pour 1L900, 3L300, et 3L900 respectivement.

## Erreur sur l'amplitude du moment dipolaire

Les erreurs sur l'amplitude du moment dipolaire pour les modèles réalistes sont représentées figure 3.7 et figure 3.10.

* Pour 1L900 l'erreur d'amplitude moyenne variait entre -50% et 0% quand la tangentialité variait entre 0° et 60°, et entre 0% et 24% quand la tangentialité variait entre 70° et 90°.

* Pour 3L300 l'erreur d'amplitude moyenne variait entre -27% et 0% quand la tangentialité variait entre 0° et 50°, et entre 0% et 6% quand la tangentialité variait entre 60° et 90°.

* Pour 3L900, l'erreur d'amplitude moyenne variait entre -11% et 0% quand la tangentialité variait entre 0° et 60°, et entre 0% et 2% quand la tangentialité variait entre 70° et 90°.

### 3.2.3 Discussion

Le but de cette simulation était d'étudier le biais introduit par différents niveaux d'approximation géométrique par rapport à un modèle réaliste de référence donné (REF). En
BRAINLOCSPH

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

1L900

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

3L300

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

3L900

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

FIG. 3.8: Erreurs de localisation pour BRAINLOCSPH et les modèles réalistes. Pour chaque position de dipôle, l'erreur de localisation correspondante est codée sous la forme d'un point coloré : rouge entre 0 et 3 mm, jaune entre 3 et 6 mm, vert entre 6 et 12 mm, et bleu au delà de 12 mm.
BRAINLOCSPH

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

1L900

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

3L300

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

3L900

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

FIG. 3.9: Erreurs d'orientation pour BRAINLOCSPH et les modèles réalistes. Pour chaque position de dipôle, l'erreur de localisation correspondante est codée sous la forme d'un point coloré : rouge entre 0 et 10°, jaune entre 10 et 20°, vert entre 20 et 40°, et bleu au delà de 40°.
<figure><img src="image_44.png" /><figcaption>FIG. 3.10: Colonne de droite : erreur de localisation (A), d'orientation (B) et d'amplitude du moment dipolaire (C) en fonction du numéro de coupe. Colonne du milieu : même chose pour les 75% de dipôles ayant les plus grandes valeurs RMS. Colonne de droite : résultats obtenus pour l'EEG par Yvert et al. [1997] pour comparaison. Dans cette étude, 3L660 était un modèle réaliste à trois couches avec 660 nœuds par couche et s'appelait IUMM. GLSPH était la sphère approchant au mieux les positions des électrodes et s'appelait Model 2. Les barres d'erreur correspondent à une erreur standard. Seuls les 1768 dipôles qui tombaient à l'intérieur de tous les maillages de cerveau et donnaient des erreurs résiduelles inférieures à 30% ont été pris en compte dans les données MEG.</figcaption></figure>
particulier, cette étude visait à établir quelles étaient les régions cérébrales pour lesquelles un
modèle sphérique donnait des résultats comparables aux modèles réalistes, et à déterminer
l'influence de certaines caractéristiques des modèles réalistes sur l'erreur de reconstruction
des dipôles. La meilleure façon d'évaluer les méthodes de localisation des sources est sans
doute l'étude de dipôles implantés. Mais on comprend aisément que pour des raisons
pratiques, peu de dipôles peuvent être considérés [Barth et al., 1986; Cohen et al., 1990;
Balish et al., 1991; Rose et al., 1991; Menninghaus et al., 1994; Leahy et al., 1998]. Les
simulations comme celle décrite dans ce chapitre nécessitent de choisir un modèle
mathématique de la tête pour créer le jeu de données, mais permettent une évaluation plus
systématique de l'erreur.

3.2.3.1 Biais dû à l'approximation sphérique

Nous avons trouvé un biais de localisation induit par l'utilisation de modèle sphérique variant de 2.5 mm en haut de la tête à 8 mm dans le lobe temporal. Cet effet de la position verticale du dipôle a également été observé en EEG sur des données simulées [Yvert et al., 1997]. Des différences de localisation importantes ont également été rapportées dans la partie inférieure de la tête entre un modèle sphérique et les plots d'une électrode implantée pour la localisation des sources de pointes ondes épileptiques [Merlet et al., 1998, Fig. 4]. Cet effet vertical est bien moins prononcé avec les modèles réalistes, que ce soit en MEG ou en EEG. Une interprétation de cet effet est que le biais de localisation induit par le modèle sphérique est plus faible dans les régions les plus sphériques de la tête qui sont également mieux couvertes par les systèmes de capteurs tête entière. L'éventail d'erreurs de localisation que nous avons trouvé est comparable à celui trouvé sur des données simulées à partir d'un modèle réaliste à une couche chez Zanow [1997] et Fuchs et al. [1998a]. Il correspond également à ce qui est trouvé chez des patients épileptiques [Balish et al., 1991; Cohen et al., 1990; Rose et al., 1991] ou sur des fantômes [Barth et al., 1986; Menninghaus et al., 1994]. A première vue, nos résultats semblent en contradiction avec ceux trouvés par Leahy et al. [1998] qui montrent qu'un modèle réaliste très fin n'améliore que de 0.5 mm l'erreur de localisation par rapport à une sphère locale. Toutefois, leurs dipôles étaient placés dans des zones imitant le sillon central et la partie supérieure du cortex occipital, or dans ces régions nous trouvons également des erreurs très similaires pour les modèles sphériques et les modèles réalistes. De plus, un de leurs dipôles (le dipôle #12) était situé dans la partie inférieure de la tête, et pour ce dipôle, le modèle réaliste utilisant la méthode des intégrales de surface donnait une erreur de 4 mm tandis que le modèle sphérique donnait une erreur d'environ 17 mm [Leahy et al., 1998, Fig. 11]. Pour les autres dipôles, situés sur la partie
supérieure de la tête, les erreurs trouvées avec le modèle réaliste et le modèle sphériques étaient similaires. Ainsi notre simulation est en accord avec leurs résultats.

Nous avons également comparé des modèles de sphères locales et globales. Pour les sphères approchant au mieux la surface du scalp, la sphère locale donnait des résultats légèrement meilleurs (1 mm) que la sphère globale dans la partie supérieure de la tête. La situation s'inversait en bas de la tête (14 mm pour LOCSPH et 7 mm pour GLSPH). En revanche, pour les sphères approchant au mieux la surface du cerveau, l'utilisation de sphères locales offrait une amélioration significative de 0.6 mm de l'erreur de localisation par rapport à une sphère globale. Les auteurs utilisant des appareils de MEG couvrant partiellement la tête (BTI 37 canaux par exemple) utilisent le plus souvent des sphères locales. Ils explorent généralement le cortex somato-sensoriel ou les aires auditives, zones pour lesquelles nous avons trouvé une légère amélioration de l'erreur avec les sphères locales approchant au mieux le scalp. Dans des études plus récentes [Mäkelä et al., 1994; Yvert et al., 1998] avec des appareils de MEG tête entière, les auteurs continuent à utiliser des sphères locales approchant au mieux le scalp et un sous-ensemble de voies voisines de la position estimée a priori des sources. Même si un tel choix se justifie pour explorer les aires auditives, il pourrait engendrer des erreurs de localisation importantes (supérieures à 8 mm) pour des sources dans le cortex frontal et inféro-temporal, à moins d'utiliser des sphères approchant au mieux la surface du cerveau. Les sphères approchant au mieux la surface du cerveau ne peuvent toutefois être utilisées que si des IRMs sont disponibles, ce qui n'est pas systématiquement le cas.

### 3.2.3.2 Influence du nombre de couches des modèles réalistes

Nous avons montré que le nombre de couches dans les modèles réaliste pourrait avoir une certaine importance pour la reconstruction des dipôles en MEG. Nous avons trouvé qu'un modèle à trois couches grossier (3L300) aboutissait à des erreurs de localisation de 2 mm inférieures à celles obtenues avec un modèle à une couche (1L900) correspondant au même temps de calcul. Ce résultat s'explique par une meilleure prise en compte des sources secondaires de courant dans le modèle réaliste à trois couches. Il ne faut toutefois pas perdre de vue que nos données simulées ont été crées avec un modèle à trois couches et non avec une description idéalement exacte des conductivités de la tête. C'est pourquoi il conviendrait de vérifier nos résultats dans des situations plus réalistes comme dans le cas de dipôles implantés chez des patients épileptiques.

### 3.2.3.3 Orientation du dipôle et amplitude du moment dipolaire

Nous avons trouvé que l'orientation du dipôle et l'amplitude du moment dipolaire pouvaient être retrouvés avec une erreur raisonnable en utilisant des modèles réalistes de la
tête. Il est généralement admis que la MEG n'est pas sensible aux composantes radiales des dipôles. Notre étude a montré que dans le cas de modèles réalistes, la distribution de champ magnétique est partiellement affectée par les composantes radiales des sources, et qu'en l'absence de bruit cette information est suffisante pour retrouver l'orientation du dipôle et l'amplitude de son moment dipolaire. Ainsi que le montre la figure 3.7 B, le modèle sphérique n'est capable de reconstruire que l'orientation et le moment dipolaire des dipôles tangentiels : l'erreur sur l'orientation et le moment dipolaire dépend donc fortement de la tangentialité. En revanche, les performances des modèles réalistes pour reconstruire l'orientation et l'amplitude du moment dipolaire ne sont que peu affectées par la tangentialité. Pour les dipôles radiaux, l'orientation est retrouvée avec une erreur inférieure à 20° pour les modèles réalistes à trois couches et 40° pour le modèle à une couche. L'amplitude du moment dipolaire est retrouvée avec une erreur inférieure à 15% pour les dipôles radiaux et 3L900. Pour les dipôles tangentiels, les modèles sphériques permettent de retrouver l'orientation des dipôles avec une erreur inférieure à 10°. Dans ce cas particulier, ils se comportent artificiellement mieux que les modèles réalistes puisque le dipôle est contraint dans ces modèles à rester tangentiel alors que l'orientation est complètement libre pour les modèles réalistes.

Il est généralement admis que les dipôles profonds et radiaux ne génèrent presque pas de champ magnétique au voisinage de la tête, or dans notre simulation comme dans d'autres données [Haueisen et al., 1995; Menninghaus et Lütkenhöner, 1995; Cuffin, 1990] de telles sources peuvent produire un champ d'amplitude égale à 20% de celui engendré par des sources tangentielles situées au même endroit. C'est pourquoi il convient d'étudier de manière plus approfondie, en particulier en présence de bruit, les capacités des modèles réalistes à retrouver les composantes radiales des sources et les sources profondes.

Par ailleurs, nous avons comparé les erreurs obtenues pour les 75% des dipôles qui correspondaient aux plus grandes valeurs RMS (figure 3.10). On constate que ce sous-ensemble de dipôles se comporte comme l'ensemble des dipôles. Ainsi, en l'absence de bruit, les erreurs de localisation et d'orientation trouvées dans cette étude dépendent davantage de la position et de la tangentialité des dipôles ainsi que du modèle géométrique que de la valeur RMS du signal.

## 3.3 Effet de l'orientation du dipôle

Nous venons de voir que les dipôles radiaux produisent des champs magnétiques dont l'amplitude n'est pas négligeable par rapport à ceux produits par des dipôles tangentiels. Nous allons ici nous intéresser plus particulièrement aux distributions de champ magnétique dans le cas de dipôles radiaux, puis nous étudierons l'effet de l'orientation sur l'amplitude des
signaux mesurés et sur les erreurs de localisation et d'orientation pour une position de dipôle fixée.

3.3.1 Méthode

Pour évaluer l'effet de l'orientation des sources, nous avons réalisé une simulation dans laquelle un dipôle était placé à une position donnée et dont l'orientation variait entre une orientation tangentielle pour le modèle sphérique associé et une orientation radiale. Pour cela, nous avons choisi une position de dipôle sur le plan médian de coordonnées (0,0,80). Nous avons alors résolu le problème direct dans un modèle réaliste pour 179 orientations du dipôle dans le plan (y-z) espacées de 1°, ce qui permettait de déduire les solutions du problème direct sur un angle de 360°. L'orientation pour laquelle le maximum sur les voies MEG de la valeur absolue du champ magnétique enregistré était minimum était considéré comme le dipôle radial. Les différentes orientations des dipôles sont représentées schématiquement figure 3.11. Elles se répartissent sur un angle de 180° avec un dipôle tous les degrés, en partant à -90° du dipôle radial pour aboutir à 90° du dipôle radial.

<figure><img src="image_45.png" /><figcaption>FIG. 3.11: Dipôles utilisés dans cette simulation : la position (0,0,80) est la même pour tous les dipôles, l'orientation reste dans le plan y-z. Les orientations se répartissent sur un angle de 180° avec un dipôle tous les degrés, en partant à -90° du dipôle dit radial pour aboutir à 90° du dipôle dit radial. Le dipôle dit radial est celui pour lequel le champ magnétique est minimal dans le modèle réaliste.</figcaption></figure>

Pour chaque orientation du dipôle, le champ magnétique a été calculé soit avec un modèle réaliste, soit avec un modèle sphérique. Le modèle sphérique a été ajusté de manière à ce que le dipôle soit radial dans le modèle sphérique quand la valeur maximale aux capteurs du champ magnétique est minimale dans le modèle réaliste, ceci afin de faciliter l'interprétation des données. Le modèle réaliste était un modèle à trois couches avec 1033, 1576, et 1459 nœuds respectivement pour le scalp, l'os et le cerveau. Le problème direct était résolu par la
méthode des intégrales de surface avec le potentiel interpolé linéairement (approche LINT) sur chaque triangle et l'approche problème isolé [Meijs et al., 1987a; Hämäläinen, 1987]. Pour 27 de ces dipôles, le problème inverse a été résolu soit avec le modèle sphérique, soit avec un modèle réaliste plus grossier composé de trois couches avec respectivement 601, 608 et 612 nœuds pour le scalp, l'os et le cerveau.

### 3.3.2 Topographie du champ magnétique engendré par un dipôle radial dans un modèle réaliste

Lors de simulations, nous avons constaté que la topographie de champ magnétique engendré par un dipôle radial n'était pas toujours constituée d'un pôle négatif et d'un pôle positif, ce qui pourrait conduire à penser que le champ magnétique en question n'est pas dû à un simple dipôle, mais plutôt à plusieurs dipôles. De tels exemples de topographies sont représentés figure 3.12. En A, la figure montre que pour le dipôle radial de notre simulation, la topographie engendrée avec le modèle réaliste contient deux pôles négatifs et deux pôles positifs, évoquant une topographie à deux dipôles. La figure 3.12 B montre cette même topographie sur une vue plus réaliste : on constate que les déformations de topographie ne sont pas simplement dûs à la projection radiale lors de la cartographie [Perrin et al., 1990]. En C et D les topographies pour deux dipôles quasi-radiaux<fnref n="1" /> sont représentées.

Dans le cas du dipôle radial dont la topographie est représentée en A et B, si l'on cherche à résoudre le problème inverse avec un modèle sphérique pour cette topographie avec un seul dipôle, on arrive à expliquer 67% du signal avec une erreur de localisation de 136 mm. En revanche, on arrive à expliquer 93 % du signal avec deux dipôles orientés à 180° l'un de l'autre, distants l'un de l'autre de 4 mm, et situés à 8 mm et 10 mm du dipôle original.

L'évolution de la topographie en fonction de l'angle par rapport au dipôle radial montre que dès que l'on se situe à 30° du dipôle radial, les topographies créées par un modèle sphérique apparaissent très différentes de celles créées par un modèle réaliste (voir figures 3.13 et 3.14). Alors que pour le modèle sphérique, la topographie est toujours identique quel que soit l'angle, et que seule l'amplitude et le signe des pôles changent, pour le modèle réaliste, on voit les pôles se déformer pour passer d'un hémisphère à l'autre et se dédoubler.

<footnote n="1">¹voir section 1.2.2 pour la définition des orientations quasi-radiales</footnote>
<figure><img src="image_46.png" /><figcaption>FIG. 3.12: A : Topographie du dipôle radial dans cette simulation. L'asymétrie de la topographie est dûe à une légère asymétrie du modèle réaliste. La position du dipôle est représentée en haut à gauche sur une vue de profil. B : Topographie du même dipôle vue en 3D. La position du dipôle est représentée en haut à gauche sur une vue de profil. C : Topographie engendrée par un dipôle quasi-radial profond situé sur un axe temporal. La position du dipôle est représentée en haut à gauche sur une vue de derrière. D : Topographie engendrée par un dipôle quasi-radial périphérique situé sur un axe temporal. La position du dipôle est représentée en haut à gauche sur une vue de derrière.</figcaption></figure>
Modèle sphérique

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

Modèle réaliste

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

FIG. 3.13: Évolution de la topographie en fonction de l'angle du dipôle : de 30° en 30°
Modèle sphérique

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

Modèle réaliste

<figure><img src="image_50.png" /><figcaption>FIG. 3.14: Évolution de la topographie en fonction de l'angle du dipôle : de 5° en 5°</figcaption></figure>
### 3.3.3 Effet de l'orientation des sources sur l'amplitude du signal mesuré

La figure 3.15 montre l'évolution du rapport entre la valeur RMS du signal et la valeur RMS pour le dipôle tangentiel en fonction de l'angle du dipôle. Pour le modèle sphérique, il est égal à la valeur absolue du sinus de l'angle. Pour le modèle réaliste, la valeur n'est pas nulle lorsque l'angle est égal à zéro, mais vaut 20% de sa valeur maximale.

En fait, si l'on pose les formules dans le repère des orientations propres du dipôle (voir section 1.2.2 pour la définition des orientations propres), on peut écrire le moment dipolaire sous la forme :

$$ \vec{M}(\theta) = \cos(\theta)\vec{M}_r + \sin(\theta)\vec{M}_t $$

où $\vec{M}_r$ et $\vec{M}_t$ sont deux orientations propres quasi-radiale et quasi-tangentielle, et $\theta$ est l'angle du dipôle. Le champ magnétique engendré est alors :

$$ \mathbf{B}(\theta) = \cos(\theta)\mathbf{B}_r + \sin(\theta)\mathbf{B}_t $$

où $\mathbf{B}_r$ et $\mathbf{B}_t$ sont les champs magnétiques engendrés par les dipôles de même position et d'orientations $\vec{M}_r$ et $\vec{M}_t$. Ils sont orthogonaux par définition. D'où :

$$ \frac{RMS(\theta)}{RMS(90)} = \sqrt{\cos(\theta)^2 \frac{RMS(B_r)^2}{RMS(B_t)^2} + \sin(\theta)^2} $$

Le rapport $\frac{RMS(B_r)}{RMS(B_t)}$ permet donc de quantifier la part de signal que le modèle sphérique ne pourra pas expliquer par rapport au modèle réaliste. Dans cette simulation, les orientations ne sont pas exactement des orientations propres puisque nous avons cherché l'orientation minimisant l'amplitude maximale du signal aux capteurs (et non le RMS) uniquement dans le plan y-z. C'est pourquoi la courbe verte de la figure 3.15 pour le modèle réaliste n'est pas exactement symétrique. Cependant, comme la tête est presque symétrique par rapport au plan y-z, ces orientations ne doivent pas être très éloignées des orientations propres. On trouve donc un rapport de 0.2 pour un dipôle à une profondeur de 38 mm. Dans une autre simulation, nous avons véritablement pris en compte les orientations propres de dipôles situés sur un axe temporal (voir figure 3.16). Le rapport variait entre 0.03 et 0.16, avec une évolution très similaire pour les deux orientations quasi-tangentielles. Ce rapport était plus élevé pour les dipôles profonds.
<figure><img src="image_51.png" /><figcaption>FIG. 3.15: Évolution du RMS en fonction de l'angle du dipôle : en rouge pour le modèle sphérique et en vert pour le modèle réaliste à trois couches. Les positions des dipôles sont représentées en vue de profil en haut à gauche.</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>FIG. 3.16: Exemple d'évolution du rapport entre le RMS du dipôle quasi-radial et les RMS des dipôles quasi-tangentiels en fonction de la position du dipôle sur un axe temporal. Les positions du dipôle sont représentées en vue de derrière en haut à gauche. La courbe bleue correspond au rapport avec le dipôle quasi-tangentiel associé à la valeur propre la plus grande et la courbe rouge au rapport avec l'autre dipôle quasi-tangentiel.</figcaption></figure>
### 3.3.4 Effet de l'orientation des sources sur les erreurs de localisation et d'orientation.

Les figures 3.17 et 3.18 montrent les erreurs de localisation et d'orientation en fonction de l'angle du dipôle. On constate que les erreurs de localisation et d'orientation se dégradent d'autant plus que le dipôle est radial pour le modèle sphérique, alors que cette dégradation est bien moindre avec le modèle réaliste.

L'erreur d'orientation très importante obtenue avec le modèle sphérique n'est pas surprenante puisque le dipôle est contraint à rester tangentiel dans ce modèle. Le modèle sphérique peut donc au mieux retrouver la composante tangentielle du dipôle. En revanche on constate une augmentation de l'erreur de localisation quand le dipôle devient radial : la composante tangentielle du dipôle ne suffit donc pas à localiser correctement la source. Les erreurs plus faibles obtenues quand l'angle du dipôle est proche de 0° doivent être mises en correspondance avec une erreur résiduelle importante pour ces angles. Comme le dipôle est initialisé avec les bons paramètres, ces erreurs faibles sont probablement dûes uniquement à la qualité de cette initialisation.

La stabilité de l'erreur obtenue avec les modèles réalistes est frappante dans cette simulation. Toutefois, cette simulation est non bruitée, et la chute de la valeur RMS quand le dipôle devient radial incite à la prudence sur la capacité des modèles réalistes à retrouver les sources radiales en présence de bruit. Cet aspect sera étudié au chapitre 4.

## 3.4 Conclusions et perspectives

La première simulation de ce chapitre montre que :

* Le biais de localisation lorsqu'on utilise un modèle sphérique croît en fonction de la position verticale de la source d'une moyenne de 2.5 mm en haut de la tête à 12 mm en bas de la tête.

* Le biais de localisation pour le modèle sphérique décroît en fonction de la tangentialité de 10 mm quand le dipôle est radial à 5 mm quand le dipôle est tangentiel.

* Il vaut mieux utiliser une sphère approchant globalement le scalp que localement si l'on ne dispose que de données de scalp et que l'on n'est pas certain que les sources se situent dans les parties supérieures ou postérieures du cerveau, et une sphère locale approchant au mieux le cerveau si l'on dispose d'IRM.

* Pour un même temps de calcul, l'utilisation d'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.
<figure><img src="image_53.png" /><figcaption>FIG. 3.17: Erreur de localisation en fonction de l'orientation du dipôle : en rouge pour le modèle sphérique et en vert pour le modèle réaliste à trois couches. En haut à droite, l'erreur résiduelle associée.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>FIG. 3.18: Erreur d'orientation en fonction de l'orientation du dipôle : en rouge pour le modèle sphérique et en vert pour le modèle réaliste à trois couches.</figcaption></figure>
* L'orientation des dipôles peut être retrouvée avec une erreur inférieure à 20° y compris pour les dipôles radiaux si on utilise un modèle réaliste.

L'étude de l'effet de l'orientation du dipôle montre que :

* Les dipôles radiaux peuvent engendrer des champs magnétiques non négligeables dont l'amplitude est égale à 20% de celle engendrée par un dipôle tangentiel à la même position et ayant la même amplitude de moment dipolaire.

* Les dipôles radiaux produisent 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 demanderaient à être confirmés par une étude en présence de bruit, ce qui sera fait au chapitre 4.

Plusieurs pistes peuvent être envisagées pour améliorer la reconstruction des sources dans les cas réalistes de données bruitées.

* Une technique de régularisation permettant d'éviter que les composantes quasi-radiales des sources soient surestimées peut être envisagée [voir Fuchs et al., 1998b].

* Le couplage de données MEG et EEG pourrait permettre de lever l'indétermination sur les composantes quasi-radiales des sources grâce aux données EEG tout en conservant la précision de localisation de la MEG, et améliorer la localisation des sources profondes.
# Chapitre 4

## Évaluation de plusieurs méthodes de couplage MEG/EEG
## 4.1 Introduction

Nous avons vu dans les sections précédentes que la MEG permettait d'obtenir des erreurs de localisation plus faibles que l'EEG dans le cas de dipôles tangentiels périphériques. Les erreurs de localisation et d'orientation en EEG sont moins dépendantes de l'orientation et de la profondeur des sources.

La comparaison des mérites relatifs des deux modalités a été le sujet d'un débat passionné [Cohen et al., 1990; Cohen et Cuffin, 1991; Hari et al., 1991; Williamson, 1991; Crease, 1991; Anogianakis et al., 1992] au terme duquel il ressort que l'EEG est plus sensible que la MEG aux variations de conductivité des milieux conducteurs, que les sources radiales et profondes sont presque silencieuses en MEG alors que l'EEG est sensible à toutes les orientations, et que ces deux modalités apportent des informations complémentaires tant pour la détection de sources neuronales que pour leur localisation. D'où l'idée d'utiliser conjointement ces deux modalités pour la résolution du problème inverse. Cette idée a été introduite bien avant [Cohen et Cuffin, 1987; Stok et al., 1987] mais n'a connu que peu de développements dans le cas de la localisation de sources dipolaires [Mosher et al., 1993; Fuchs et al., 1998b]. Des techniques de couplage adaptées aux modèles de sources distribuées ont été étudiées [Hasson, 1991; Baillet et al., 1999], mais elles font intervenir des techniques très différentes puisqu'il s'agit dans ce cas de résoudre un système linéaire sous-déterminé. La possibilité offerte par le logiciel CURRY® de réaliser le couplage avec des modèles réalistes a permis un début de validation expérimentale du couplage dans le cas de données évoquées somesthésiques pour de simples sources dipolaires [Fuchs et al., 1998b] ou pour des sources distribuées [Babiloni et al., 1999a;b; Kasai et al., 1999]. Une utilisation expérimentale d'une technique de couplage avec un modèle sphérique (logiciels BESA® et Neuromag®) et des données évoquées auditives est également présentée dans Huotilainen et al. [1998].

Plusieurs méthodes de couplage sont proposées dans la littérature : une partie d'entre elles utilisent séquentiellement les données MEG et EEG. Ainsi, Cohen et Cuffin [1987] proposent une méthode en trois étapes :

1. Le problème inverse MEG est résolu par la méthode habituelle, il donne pour solution une configuration de sources $d_{MEG}$. Ces sources devraient être tangentielles et superficielles.

2. Le problème direct EEG est résolu pour la configuration de sources $d_{MEG}$ et donne un jeu de valeurs de potentiel aux électrodes $V(d_{MEG})$ qui sera soustrait aux valeurs mesurées d'EEG ($V$).

3. Le problème inverse EEG est résolu pour $V - V(d_{MEG})$ et donne pour solution une configuration de sources $d_{ EEG}$. Ces sources devraient être soit profondes, soit radiales.