\documentclass[XUPS,XML,SOM,Unicode,francais, NoFloatCountersInSection, NoEqCountersInSection]{cedram}
\usepackage{xups22-03}
\setcounter{tocdepth}{2}
\mathtoolsset{showonlyrefs}
%\XUPScorrections

\hyphenation{décrits magné-ti-que magné-ti-ques besoin ortho-go-nal ortho-go-na-le fini finis}

\begin{document}
\frontmatter
\title[Modélisation du signal d'IRM de diffusion par une EDP]{Modélisation du signal d'IRM de\nobreakspace diffusion par une\nobreakspace équation aux\nobreakspace dérivées partielles}

\author[\initial{J.-R.} \lastname{Li}]{\firstname{Jing-Rebecca} \lastname{Li}}
\address{Equipe IDEFIX, INRIA Saclay,
Unité de Mathématiques Appliquées (UMA), ENSTA Paris,
828, Boulevard des Maréchaux, 91762 Palaiseau, France}
\email{jingrebecca.li@inria.fr}
\urladdr{http://www.cmap.polytechnique.fr/~jingrebeccali/}

\thanks{Journées X-UPS 2022. Mathématiques et biologie. Éditions de l'École polytechnique, 2022}

\begin{abstract}
Le signal d'IRM (pondéré en diffusion) est la somme de la magnétisation dans un volume de tissu cellulaire (un voxel). La magnétisation à l’échelle de la micro-structure cellulaire peut être modélisée par une équation aux dérivées partielles qui s'appelle équation de Bloch-Torrey. Ce qui rend difficile la solution numérique de cette équation est la présence d'interfaces complexes (à savoir, les membranes des cellules) sur lesquelles la solution est discontinue. Dans ce texte on va découvrir le comportement de la solution dans des géométries qui représentent les cellules cérébrales ainsi que ce que les scientifiques espèrent pouvoir apprendre grâce à cette modalité d'imagerie.
\end{abstract}

\maketitle
\vspace*{-1.5\baselineskip}
\tableofcontents

\mainmatter
\Changel

\vspace*{-2.5\baselineskip}\vskip0pt
\section{Introduction}

L'imagerie par résonance magnétique (IRM) repose sur le phénomène de résonance magnétique nucléaire (RMN) observé indépendamment par F.\,Bloch \cite{Bloch46} et E.\,Purcell \emph{et al.} \cite{P-T-P46} en 1946. Le principe de base peut être décrit succinctement de la façon suivante. Un champ magnétique principal constant de forte intensité est appliqué au tissu biologique étudié dans une direction fixée (parallèle au scanner d'IRM). Suite à l'application d'impulsions de champs magnétiques transverses au champ magnétique principal, les spins des atomes d'hydrogène contenus dans les molécules d'eau changent d'orientation, puis reviennent à leur position initiale en émettant un signal électrique. La vitesse avec laquelle ces atomes reviennent à leur orientation naturelle est fonction de différentes propriétés des organes étudiés. Des détecteurs analysent les infimes champs électriques créés par ces atomes d’hydrogène et les utilisent pour recréer une image en 2D (une \og tranche\fg) ou en 3D (en compilant les tranches). Les premières reconstructions d'images d'objets à partir du signal d'IRM sont réalisées en 1973 par P.\,Lauterbur. En 1975, R.\,R.\,Ernst propose d'utiliser la transformation de Fourier pour analyser le codage du signal IRM et la même année, P.\,Mansfield produit les premières images de tissus humains. En 1977, R.\,V.\,Damadian réalise les premières images IRM d'un corps humain vivant. Le lecteur trouvera dans l'article \cite{Goldman04} plus de détails sur cette histoire.\enlargethispage{2\baselineskip}%

Dans ces notes, nous nous intéressons à un procédé d'IRM plus sophis\-tiqué, appelé l'IRM pondéré en diffusion, qui permet l'étude des mouvements microscopiques des molécules d'eau (mouvement brownien qui caractérise la diffusion de l'eau). Principalement utilisée pour l’imagerie cérébrale, l'IRM de diffusion est l’unique méthode de cartographie des connexions neuronales in vivo. Cette technique de diffusion consiste à modifier le champ magnétique principal en un champ magnétique non uniforme (c'est-à-dire dépendant de la~position), en décalage dans le temps avec les champs magnétiques d'imagerie \hbox{décrits} précédemment.

La magnétisation, aussi appelée aimantation, est la somme macro\-scopique des (moments magnétiques de) spins des atomes d'hydrogène de l'eau. Elle est représentée par un vecteur qui dépend du temps et de la position. À l'échelle de la micro-structure cellulaire, la \emph{magné\-tisation transverse induite par le champ non uniforme (projection orthogonale de la magnétisation dans le plan transverse au champ magnétique principal)} peut être modélisée par une équation aux déri\-vées partielles, appelée l'équation de Bloch-Torrey \cite{Torrey56}. Ce~qui rend difficile la résolution numérique de cette équation est la présence d'interfaces complexes (c'est-à-dire les membranes des cellules) sur lesquelles la solution est discontinue.

Nous présenterons le modèle mathématique permettant de décrire la magnétisation transverse et le signal d'IRM qui en découle, et nous étudierons des approximations numériques de ceux-ci. Nous découvrirons le comportement de la solution dans des géométries qui représentent les cellules cérébrales et nous expliquerons ce que les scientifiques espèrent pouvoir apprendre grâce à cette modalité d'imagerie.\enlargethispage{\baselineskip}%

\section{Modèle mathématique}

Nous considérons un petit volume de tissu biologique cérébral, défini comme un domaine tridimensionnel $\Omega = \bigcup_{i = 1}^{N_\textup{cmpt}} \Omega_i $ composé de $N_\textup{cmpt}$ compartiments $\{\Omega_i\}_{1 \leq i \leq N_\textup{cmpt}}$, où les compartiments sont des parties fermées, d'intérieurs deux à deux disjoints qui correspondent à des cellules, à des parties de cellules, ou à un espace extra-cellulaire. L'interface perméable entre deux compartiments est notée $\Gamma_{i j} = \Omega_i \cap \Omega_j$ pour $i \neq j$, $(i, j) \in \{1, \dots, N_\textup{cmpt}\}^2$. Enfin, soit $\partial \Omega$ le bord extérieur du domaine, et $\Gamma_i = \Omega_i \cap \partial \Omega$ sa restriction à $\Omega_i$.

Si le volume est constitué de $N_\textup{cell}$ cellules entourées par un espace extra-cellulaire (ECS), l'ECS sera le dernier compartiment. Pour les cellules sphériques, seul l'ECS a un bord extérieur, et toutes les cellules ont une interface avec l'ECS. Les cellules peuvent également contenir des noyaux, et dans ce cas, plusieurs compartiments représenteront la cellule. Pour les cellules cylindriques, le haut et le bas de chaque cylindre constituent un bord extérieur. Les cellules peuvent avoir des couches de myéline qui influent sur la perméabilité entre elles.

\subsection*{Équation de Bloch-Torrey}
Nous précisons tout d'abord, de manière simplifiée, les différents objets issus de la physique intervenant dans l'EDP de Bloch-Torrey. On rappelle que le spin est une propriété quantique d'une particule qui peut être caractérisée de différentes \hbox{façons}. Ici, nous caractérisons le spin du proton d'un atome d'hydro\-gène par son moment magnétique $\mu\in\R^3$, et nous confondrons ces deux notions, comme ceci est fréquemment l'usage. La magnétisation, ou aimantation, est la somme macroscopique des spins des protons des atomes d'hydrogène des molécules d'eau contenues dans le volume considéré.

Le volume $\Omega$ est tout d'abord soumis à un champ magnétique principal $\mathcal B_0$ constant de forte intensité $B_0$ dans une seule direction, que nous noterons $Oz$. Soumis à ce champ magnétique, les spins des atomes d'hydrogène s'orientent parallèlement à $\mathcal B_0$ et suivent un mouvement de rotation (appelé précession) autour de l'axe $0z$ avec une fréquence $\omega_0=\gamma B_0$ (appelée fréquence de Larmor), où $\gamma = 2,67513 \times 10^8 \, \mathrm{rad \, s^{-1} T^{-1}}$ est le rapport gyromagnétique du proton de l'eau.\enlargethispage{-2\baselineskip}%

Le volume $\Omega$ est ensuite soumis à des impulsions de champs magné\-tiques dans des directions transverses à $Oz$, c'est-à-dire dans le plan $Oxy$.

En IRM de diffusion, on applique, en supplément aux champs décrits précédemment, un gradient de champ magnétique au tissu pour coder la diffusion de l'eau. Ceci signifie que l'on modifie, selon des modalités dans le temps bien déterminées, le champ magnétique principal $\mathcal B_0$ en un champ magnétique non uniforme, dit de diffusion, en ajoutant un champ magnétique de la forme ${\mathcal B}_{\vec{g}}=(0,0,\vec{g}\cdot\vec{x})$ où $\vec{g}\in\R^3$ est appelé le gradient du champ magnétique de diffusion, ou plus simplement, gradient de diffusion.

La façon dont on applique le gradient de champ magnétique de diffusion au cours du temps est codé par une fonction $f$ du temps, appelée profil temporel du gradient de diffusion, ou encore séquence de codage de diffusion.

Au cours des différentes séquences d'applications de champs ma\-gnétiques décrites précédemment, la magnétisation $\mathcal M_{i}$ correspondant au compartiment $\Omega_i$ du tissu biologique étudié évolue en fonction de la position $\vec{x}$ et du temps $t$.

En notant $(\mathcal M_{i,x}, \mathcal M_{i,y})$ ses coordonnées dans un repère fixe $Oxy$ du plan transverse à $Oz$, on définit
la fonction $$M_i(\vec{x},t)=(\mathcal M_{i,x}(\vec{x},t)+\imunit \mathcal M_{i,y}(\vec{x},t))e^{\imunit \gamma B_0 t}$$ où $\imunit$ est l'unité imaginaire ($\imunit^2=-1$). Dans toute la suite de ce texte, nous appellerons cette fonction $M_i$ la \emph{magnétisation transverse}, ou tout simplement la \emph{magnétisation}, du compartiment $\Omega_i$. Cette \hbox{fonction} mesure l'impact de l'ajout du champ magnétique $\mathcal B_{\vec{g}}$ sur la composante de la magnétisation dans le plan transverse.
Elle dépend du gradient de diffusion $\vec{g}$ et du profil temporel $f$ et elle satisfait à l'EDP de Bloch-Torrey \cite{Torrey56}:
\begin{multline}\label{EDPBT}
\frac{\partial}{\partial t} M_i(\vec{x}, t) = -\imunit \gamma f(t) \vec{g} \cdot \vec{x} \, M_i(\vec{x}, t) + \nabla \cdot D_i \nabla M_i(\vec{x}, t), \\
\vec{x} \in \Omega_i, \ i \in \{1, \dots, N_\textup{cmpt}\},
\end{multline}
où $\gamma$ est comme plus haut et $D_i$ est le coefficient de diffusion intrinsèque dans le compartiment $\Omega_i$.

Le deuxième terme du membre de droite est spécifique au procédé de diffusion et provient de l'étude des forces exercées par le gradient de diffusion sur les spins.

On désignera par $M(\vec{x}, t)$ la \emph{magnétisation globale} définie pour presque tout $\vec{x} \in \Omega$ par $M(\vec{x}, t)=M_i(\vec{x},t)$ lorsque $\vec{x}$ est un point intérieur de $\Omega_i$. Notons que la magnétisation peut être discontinue aux interfaces $\Gamma_{i j}=\Omega_i\cap \Omega_j$.

Outre l'EDP \eqref{EDPBT}, l'équation de Bloch-Torrey comporte des conditions initiales et des conditions sur les bords des compartiments, appe\-lées également conditions aux limites.

Les conditions initiales sont supposées constantes à l'intérieur de chaque compartiment:
\begin{equation}
M_i(\vec{x}, 0) = \rho_i, \quad \vec{x} \in \mathrm{int}(\Omega_i), \quad i \in \{1, \dots, N_\text{cmpt}\},
\end{equation}
avec $\rho_i\geq 0$. (Il est toujours possible de choisir un repère $Oxy$ du plan transverse à $Oz$ de sorte que la magnétisation au temps $t=0$ vérifie cette propriété.)

Les conditions aux limites extérieures sont données par
\begin{equation} \label{eq:outer_boundary}
D_i \nabla M_i(\vec{x}, t) \cdot \vec{n}_i(\vec{x}) = 0, \quad \vec{x} \in \Gamma_i=\Omega_i\cap \partial\Omega,\ i \in \{1, \dots, N_\text{cmpt}\},
\end{equation}
où $\vec{n}_i$ est le vecteur normal unitaire pointant vers l'extérieur du compartiment $\Omega_i$. Ces conditions traduisent le fait que l'eau ne se diffuse pas à l'extérieur de $\Omega$. Notez que nous pouvons avoir $\Gamma_i = \emptyset$, car tous les compartiments ne touchent pas nécessairement le bord extérieur.

L'équation de Bloch-Torrey doit également être complétée par des conditions d'interface. Nous rappelons que l'interface entre $\Omega_i$ et $\Omega_j$ est $\Gamma_{i j}$. La première condition d'interface traduit la continuité du flux de la magnétisation sur le bord $\Gamma_{ij}$ (le flux des spins de $\Omega_i$ vers $\Omega_j$ est l'opposé du flux des spins de $\Omega_j$ vers $\Omega_i$). La deuxième condition exprime le comportement du flux de la magnétisation selon la nature de la membrane $\Gamma_{ij}$. Elle dépend du coefficient de perméabilité \hbox{$\kappa_{i j}\!\geq\!0$} à travers $\Gamma_{i j}$:
\begin{align*}
D_i \nabla M_i(\vec{x}, t) \cdot \vec{n}_i(\vec{x}) & = -D_j \nabla M_j(\vec{x}, t) \cdot \vec{n}_j(\vec{x}), \\
&\hspace*{1.5cm}\notag\vec{x} \in \Gamma_{i j}, \ (i, j) \in \{1, \dots, N_\textup{cmpt}\}^2, & \\
D_i \nabla M_i(\vec{x}, t) \cdot \vec{n}_i(\vec{x}) & = \kappa_{i j} \left(c_{i j} M_j(\vec{x}, t) - c_{j i} M_i(\vec{x}, t)\right), \\
&\hspace*{1.5cm}\notag\vec{x} \in \Gamma_{i j}, \ (i, j) \in \{1, \dots, N_\textup{cmpt}\}^2.
\end{align*}

Les conditions aux limites ci-dessus sont appelées conditions de Neumann car elles portent sur le gradient de la magnétisation.

Ici, le coefficient de perméabilité, dont on suppose qu'il est positif ou nul, est une caractéristique de la membrane uniquement (\hbox{$\kappa_{i j} \!=\! \kappa_{j i}$}). Si $\kappa_{i j} = 0$, l'interface $\Gamma_{i j}$ n'est pas perméable, et la condition est similaire à celle au bord extérieur dans l'équation \eqref{eq:outer_boundary}.

Les deux poids $c_{i j}$ et $c_{j i}$ tiennent compte de la répartition des atomes d’hydrogène entre les deux compartiments. Ils peuvent tous deux être fixés à $1$, auquel cas une répartition uniforme entre les compartiments est favorisée en l'absence de gradient. Pour des répartitions initiales différentes, ce qui se traduit par \soutnew{$\rho_i$ et $\rho_j$}{$\rho_i \neq \rho_j$}, nous autorisons également des poids non symétriques, par exemple $c_{i j} = \spfrac{2 \rho_i}{\rho_i + \rho_j}$ et $c_{j i} = \spfrac{2 \rho_j}{\rho_i + \rho_j}$ comme proposé dans \cite{Lee2021}. Cela garantit que la répartition initiale non uniforme est préservée si le gradient $\vec{g}$ est nul. Le coefficient de normalisation $2 / (\rho_i + \rho_j)$ garantit que $c_{i j} = c_{j i} = 1$ si $\rho_i = \rho_j$.

Un profil temporel couramment utilisé est la séquence d'écho de spin à gradient pulsé (PGSE) \cite{Stejskal1965}, avec deux impulsions rectangulaires de durée $\delta$, séparées par un intervalle de temps $\Delta - \delta$, pour lequel le profil $f$ est donné pour $ t\in[0,T_\textup{e}]$ par
\begin{equation*}
f(t) =
\begin{cases}
1, \quad & t_1 \leq t \leq t_1 + \delta, \\
-1, \quad & t_1 + \Delta < t \leq t_1 + \Delta + \delta, \\
0, \quad &\text{sinon},
\end{cases}
\end{equation*}
où $t_1$ est le temps de départ de la première impulsion de gradient et~$T_\textup{e}$, appelé le temps d'écho, est le temps auquel le signal est mesuré (avec $t_1 + \Delta \geq T_\textup{e} / 2$) (voir la figure \ref{fig:PGSE}). Dans ce qui suit, nous fixerons $t_1 = 0$ et $T_\textup{e} = \Delta + \delta$ pour faire simple.
\begin{figure}[!ht]
\centering
{\includegraphics[width=0.6\columnwidth]{PGSEschema.png}}
\caption{Schéma séquentiel de l'écho de spin à gradient pulsé (PGSE). {\bf RF}: le profil temporel de l'application des impulsions de radio-fréquence.
$\mathbf{G(t)}$: profil temporel de l'appli\-cation des gradients de codage de diffusion.
$\mathbf{f(t)}$: profil temporel effectif de l'application des gradients de codage de la diffusion,
en tenant compte de l'impulsion de 180 degrés.
$\mathbf {\delta}$: durée des impulsions du gradient de codage de diffusion. $\mathbf{\Delta}$: délai entre le début de $2$ impulsions consécutives d'une séquence temporelle.
{\bf Écho}: moment où le signal est acquis.}
\label{fig:PGSE}
\end{figure}

Le signal IRM de diffusion est mesuré au temps d'écho $t = T_\textup{e}$. Ce~signal est l'intégrale spatiale de la magnétisation finale $M(\vec{\cdot}, T_\textup{e})$:
\begin{equation*}
S = \int_{\vec{x} \in \Omega} M(\vec{x}, T_\textup{e}) \, \mathrm{d}\bx\quad(\text{$\mathrm{d}\bx$ est la mesure de Lebesgue usuelle}).
\end{equation*}

Dans une expérience d'IRM de diffusion, la séquence d'impulsions (profil temporel $f$) est généralement fixe, tandis que l'on fait varier l'amplitude (et éventuellement la direction) de $\vec{g}$. Le signal $S$ est généralement représenté par rapport à une quantité appelée la valeur de $b$. La valeur de $b$, qui dépend de $\vec{g}$ et de $f$, est définie comme suit
\begin{equation*}
b(\vec{g}, f) = \gamma^2 \|\vec{g}\|^2 \int_0^{T_\textup{e}} \soutnew{\, \mathrm{d} s}{} \left(\int_0^s f(t) \, \mathrm{d} t\right)^2 \soutnew{}{\, \mathrm{d} s}.
\end{equation*}
Pour le PGSE, la valeur de $b$ est donnée par \cite{Stejskal1965}:
\begin{equation*}
b(\vec{g}, \delta, \Delta) = \gamma^2 \|\vec{g}\|^2 \delta^2(\Delta - \delta / 3).
\end{equation*}
Ces définitions s'expliquent par le fait que dans un milieu homogène (c'est-à-dire $\Omega=\mathbb R^3$), l'atténuation du signal est de $\euler^{-D b}$, où $D$ est le coefficient de diffusion intrinsèque.

\section{Modèle géométrique}

Notre objectif est de trouver d'une part, une solution approchée de l'équation de Bloch-Torrey et, d'autre part, d'étudier plus précisément le signal obtenu. Pour cela, nous utilisons une méthode de discrétisation spatiale, que nous décrirons à la section \ref{sect-elts finis}. Cette méthode est basée sur un maillage tétraédrique du tissu biologique étudié, qui est obtenu de la façon suivante.

\begin{figure}[htb]
	\centering
	\includegraphics[width=0.4\textwidth]{02a_pyramidal2aFI_medical_seg_zoom} \quad
	\includegraphics[width=0.4\textwidth]{02a_pyramidal2aFI_wrapped_mesh_zoom}
	\caption{\emph{Gauche}: une description de surface d'un neurone pyramidal publiée dans \emph{NeuroMorpho.Org} \cite{ascoli9247}
qui contient
de nombreuses auto-intersections et proximités; elle ne peut pas être utilisé
pour la génération du maillage des éléments finis.
\emph{Droite}: une nouvelle triangulation de surface qui répare les auto-intersections et les proximités; elle est prête à être utilisée pour les opérations suivantes de
génération de maillage par éléments finis.}
	\label{fig:02a_pyramidal2aFI_image}
\end{figure}

Nous considérons les reconstructions morphologiques (fichiers au format SWC) publiées dans \emph{NeuroMorpho.Org} \cite{ascoli9247},
la plus grande collection de reconstructions neuronales en 3D accessibles au public.
Ces descriptions de surface des neurones, obtenues à partir d'images captées par des caméras à très haute précision, ne sont pas d'une qualité suffisante pour des calculs numériques
car elles contiennent de nombreuses auto-intersections et proximités (voir figure \ref{fig:02a_pyramidal2aFI_image}, à gauche).

Nous avons utilisé un logiciel \cite{ansa} pour corriger et améliorer la qualité des descriptions de la surface des neurones.
Nous avons ensuite produit de nouvelles triangulations de surface (voir figure \ref{fig:02a_pyramidal2aFI_image}, à droite) qui sont prêtes à être utilisées pour les calculs numériques.
Les nouvelles triangulations de surface sont traitées par un autre logiciel (GMSH \cite{Geuzaine2009}) pour obtenir des
mailles tétraédriques de volume.\enlargethispage{\baselineskip}

\section{Solution numérique par éléments finis}\label{sect-elts finis}

La discrétisation spatiale de l'équation de Bloch-Torrey est basée sur la méthode des éléments finis \cite{Nguyen2014}, dont nous allons décrire les principales étapes ci-dessous. L'approximation de la magnétisation satisfait à une équation différentielle ordinaire (EDO) par rapport au temps, en dimension finie.
La résolution dans le temps de cette équation s'obtient ensuite par une méthode de résolution d'EDO de type Runge-Kutta.

Le modèle géométrique décrit dans la section précédente consiste donc à approcher le tissu biologique $\Omega$ et les compartiments $\Omega_i$ par une union finie de tétraèdres. Le principe général de la méthode des éléments finis est d'approcher la solution de l'équation considérée par une fonction d'un espace vectoriel de dimension finie formé de fonctions polynomiales par morceaux, appelé espace des éléments finis.

L'espace dans lequel nous chercherons une approximation de la solution de l'EDP est étroitement lié à la discrétisation de $\Omega$. Pour simplifier la présentation des différents résultats obtenus dans notre situation, nous supposons dans toute la suite de ce texte que les membranes entre les différents compartiments intérieurs sont imperméables (c'est-à-dire $\kappa_{i,j}=0$ pour toute interface $\Gamma_{ij}$). Les simulations numériques présentées dans les deux dernières parties concernent des neurones qui satisfont à cette propriété.

\subsection{L'espace des éléments finis}
Sous l'hypothèse simplificatrice d'imperméabilité des membranes intérieures, nous n'avons plus besoin de considérer tous les compartiments en même temps, nous pouvons considérer un compartiment à la fois, car il n'interagit avec aucun autre. La condition aux limites sur le compartiment $\Omega_i$ devient simplement
\[
D_i \nabla M_i(\vec{x}, t) \cdot \vec{n}_i(\vec{x}) = 0, \quad \vec{x} \in \partial \Omega_i, \ i\in \{1, \dots, N_\textup{cmpt}\}.
\]

Le modèle géométrique a permis d'approcher le domaine $\Omega_i$ par l'union de $N_T$ éléments tétraédriques:
\begin{equation*}
\Omega_i \approx \calT^h \equiv \bigcup_{l = 1}^{N_T} T_l,
\end{equation*}
où $T_l$ est le $l$-ème tétraèdre, appelé élément fini. L'union $\calT^h$ des éléments finis constitue le maillage par éléments finis du domaine~$\Omega_i$. Le nombre $h$ indique la taille des éléments. Il est donc possible d'avoir plusieurs maillages de tailles différentes pour $\Omega_i$, où $h$ et $N_T$ sont différents. Plus la taille des éléments (finis) est petite, plus la solution approximative est précise.

L'espace de fonctions associé à un maillage $\calT^h$ que nous utiliserons est un sous-espace de l'espace fonctionnel le plus simple, noté~$\bbP_1$, formé des fonctions affine sur chaque~tétraèdre, globalement continues.

Pour définir une base de cet espace, nous introduisons des points particuliers, appelés n{\oe}uds, ou points par éléments finis du mail\-lage~$\calT^h$. Ici, les nœuds sont les sommets de tous les éléments \hbox{finis} (c'est-à-dire des tétraèdres), soit $4$ n{\oe}uds pour chaque élément. Cepen\-dant, il est important de noter que le nombre de nœuds est bien inférieur à $4 N_T$ car les éléments se touchent, ainsi le même nœud peut appartenir à plusieurs éléments.

On note $v_l^1,v_l^2,v_l^3,v_l^4$ les sommets de l'élément $T_l$. Ainsi, l'ensemble des nœuds est
\begin{equation*}
\{P_1,P_2,\dots, P_{N_P} \}= \bigcup_{l=1}^{N_T} \{v_l^1,\dots, v_l^4\}.
\end{equation*}
où $N_P$ est le nombre total de n{\oe}uds dans $\calT^h$.

On définit un ensemble de fonctions de base
\begin{equation*}
\varphi_j(x,y,z),\quad j = 1,\dots, N_P,
\end{equation*}
de la manière suivante. La fonction de base $\varphi_j(x,y,z)$ est la fonction affine par morceaux en $x$, $y$, $z$ définie par:
\begin{multline}\label{eq:bf_support}
\varphi_j(x,y,z)\\
= \begin{cases}
0 & \text{sur } T_l \text{ si } P_j \not\in \{v_l^1,\dots, v_l^4\}\\
&(P_j \text{ n'est pas un n{\oe}ud de } T_l),\\
a_j^lx+b_j^ly+c_j^lz+d_j^l & \text{sur } T_l \text{ si } P_j \in \{v_l^1,\dots, v_l^4\}\\
&(P_j \text{ est un n{\oe}ud de } T_l).
\end{cases}
\end{multline}
dont les coefficients $a_j^l$, $b_j^l$, $c_j^l$, $d_j^l$ sont déterminés par les 4 contraintes suivantes imposées aux sommets de $T_l$:
\begin{equation}\label{eq:bf_constraints}
\varphi_j(x,y,z) = \begin{cases} 1 & \text{si } (x,y,z) = P_j \\
0 & \text{si } (x,y,z) \neq P_j \text{ et } (x,y,z) \in \{v_l^1,\dots, v_l^4\}.
\end{cases}
\end{equation}
La fonction $\varphi_j$ est donc uniquement déterminée par le n{\oe}ud $P_j$.

En bref, pour obtenir la fonction de base $\varphi_j(x,y,z)$, il faut trouver tous les $T_l$ pour lesquels $P_j$ est un sommet et trouver les coefficients $a_j^l$, $b_j^l$, $c_j^l$, $d_j^l$ qui définissent $\varphi_j$ sur l'élément $T_l$ en imposant les contraintes de l'équation \eqref{eq:bf_constraints}. Le support de $\varphi_j$ est petit, comme on peut le voir dans l'équation \eqref{eq:bf_support}, $\varphi_j$ est non nulle uniquement sur les éléments tétraédriques dont $P_j$ est un sommet. Sur la grande majorité des éléments, $\varphi_j$ est identiquement nulle.

\subsection{Formulation faible de l'EDP}
\label{sec:fem}

Comme nous supposons que les interfaces sont imperméables, nous pouvons considérer l'équation dans un seul compartiment $\Omega_i$ à la fois, et donc l'équation étudiée~est:
\begin{equation*}%\label{EDPBT_onecmpt}
\begin{aligned}
&\frac{\partial}{\partial t} M_i(\vec{x}, t) \!=\! -\imunit \gamma f(t) \vec{g} \cdot \vec{x} \, M_i(\vec{x}, t) \!+\! \nabla \cdot D_i \nabla M_i(\vec{x}, t), \ \vec{x} \in \Omega_i, \\
&D_i \nabla M_i(\vec{x}, t) \cdot \vec{n}_i(\vec{x}) = 0, \quad \vec{x} \in \partial\Omega_i.
\end{aligned}
\end{equation*}

La méthode des éléments finis ne résout pas la formulation forte ci-dessus, mais plutôt une formulation faible dans l'espace de Sobolev $H^1(\Omega_i)$ de l'EDP, où n'interviennent que des dérivées spatiales du premier ordre (et pas du second ordre) de la solution $u$. On rappelle que $H^1(\Omega_i)$ est l'espace des fonctions $v\in L^2(\Omega_i)$ dont les dérivées partielles d'ordre $1$ (au sens des distributions) appartiennent à $L^2(\Omega_i)$.

La formulation faible est obtenue comme suit: nous multiplions la formulation forte par une fonction test $v$, c'est-à-dire une fonction quelconque de $H^1(\Omega_i)$, nous intégrons sur le domaine $\Omega_i$ et nous utilisons l'identité de Green. Plus précisément,
en notant $\mathrm{d}\bx$ la mesure volumique du domaine $\Omega_i$, c'est-à-dire $\mathrm{d}\bx = dx \,dy \,dz$, on obtient
\begin{multline*}%\label{EDPBT_onecmpt_weak}
\int_{\Omega_i} \frac{\partial}{\partial t} M_i(\vec{x}, t) v\, \mathrm{d}\bx= \int_{\Omega_i}-\imunit \gamma f(t) \vec{g} \cdot \vec{x} \, M_i(\vec{x}, t) v\, \mathrm{d}\bx\\
+ \int_{\Omega_i}\nabla \cdot D_i \nabla M_i(\vec{x}, t)v\, \mathrm{d}\bx.
\end{multline*}
L'identité de Green donne:
\begin{multline*}%\label{eq:green}
\int_{\Omega_i} -\lap{D_i}{M_i}v\; d \bx = \int_{\Omega_i} \left (D_i\nabla M_i\right) \cdot \nabla v\, \mathrm{d}\bx\\
- \int_{\partial \Omega_i} \left (D_i\nabla M_i \right) \cdot \bn_i \; v\; \mathrm{d}\bs,
\end{multline*}
où $d\bs$ indique l'intégration de surface, sur la limite du domaine $\Gamma_i = \partial \Omega_i$.
Comme la perméabilité est nulle, l'intégrale de surface ci-dessus est nulle et nous obtenons l'équation qui ne porte que sur $M_i$ et $\nabla M_i$:
\begin{multline}\label{eq:weak4}
\frac{\partial}{\partial t} \int_{\Omega_i} M_i\,v \; \mathrm{d}\bx+\int_{\Omega_i} \left (D_i\nabla M_i\right) \cdot \nabla v\,\mathrm{d} \bx\\
+\int_{\Omega_i} \imunit \gamma f(t) \vec{g} \cdot \vec{x} \, \,M_i\,v \; \mathrm{d}\bx = 0.
\end{multline}

\subsection{Discrétisation par éléments finis}

L'espace de fonctions dans lequel nous cherchons une solution approchée de la magnétisation restreinte au compartiment $\Omega_i$ est l'espace des fonctions affines par morceaux sur les éléments tétraédriques, globalement continues, espace défini par:
\begin{equation*}\textstyle
\calU^h = \left \{u(x,y,z,t) = \sum_{j=1}^{N_P} U_j(t) \; \varphi_j(x,y,z), \ U_j(t) \in \R \right \}.
\end{equation*}
L'exposant $h$ rappelle le maillage d'éléments finis sous-jacent $\calT^h$ sur lequel les fonctions de base sont définies.

Par construction des $\varphi_j$, le coefficient $U_j$ d'une fonction $u\in\calU^h$ est juste la valeur de $u$ sur le nœud $P_j$:
\begin{equation*}
u(P_j,t) = U_j(t).
\end{equation*}
Il s'agit d'une propriété très utile et pratique du choix de l'espace de fonctions des éléments finis que nous avons choisi.

Nous cherchons donc une solution approchée de la forme
\begin{equation*}%\label{eq:u}
u^h(x,y,z,t) = \sum_{j=1}^{N_p} U_j(t)\ \varphi_j(x,y,z), \quad U_j(t) \in \R.
\end{equation*}

Pour trouver la solution approchée $u^h$, il suffit de trouver les coef\-ficients $U_j(t)$. Et puisque $U_j$ coïncide avec la valeur de la fonction au nœud $P_j$, nous obtenons en même temps les valeurs de la solution approchée aux nœuds du maillage des éléments finis.

Nous appliquons l'équation \eqref{eq:weak4} à la solution approchée $u^h$ et, pour une fonction test $v\in \calU^h$, nous obtenons:
\begin{equation*}%\label{eq:weak5a}
\frac{\partial}{\partial t} \int_{\Omega_i} u^h\,v \; \mathrm{d}\bx+\int_{\Omega_i} \left (D_i\nabla u^h\right) \cdot \nabla v\,\mathrm{d} \bx
+\int_{\Omega_i} \imunit \gamma f(t) \vec{g} \cdot \vec{x} \, \,u^h\,v \; \mathrm{d}\bx = 0.
\end{equation*}
et donc $u^h$ satisfait à l'équation de contrainte
\begin{multline*}
\sum_{j=1}^{N_P} \frac{\partial}{\partial t}U_j(t) \;\int_{\Omega_i} \, \varphi_j\,v \; \mathrm{d}\bx
+\sum_{j=1}^{N_P} U_j(t) \int_{\Omega_i}\left (D_i \nabla \; \varphi_j\right) \cdot \nabla v\,\mathrm{d} \bx\\
+ \sum_{j=1}^{N_P} U_j(t) \;\int_{\Omega_i} \imunit \gamma f(t) \vec{g} \cdot \vec{x} \, \, \varphi_j\,v \; \mathrm{d}\bx
= 0.
\end{multline*}
Comme il y a $N_P$ coefficients inconnus (les $U_j(t)$), nous avons besoin de $N_P$ équations de contrainte. Pour cela, nous prenons comme choix de fonction test $v\in \calU^h$ chacune des fonctions de base $\varphi_i$, \hbox{$i = 1,\dots,{N_P}$}. Les $N_P$ équations de contrainte sont:
\begin{multline*}
\sum_{j=1}^{N_P} \frac{\partial}{\partial t}U_j(t) \;\int_{\Omega_i} \, \varphi_j\,\varphi_i \; \mathrm{d}\bx
+\sum_{j=1}^{N_P} U_j(t) \int_{\Omega_i} \left (D_i \nabla \; \varphi_j\right) \cdot \nabla \varphi_i\,\mathrm{d} \bx\\
+ \sum_{j=1}^{N_P} U_j(t) \;\int_{\Omega_i} \imunit \gamma f(t) \vec{g} \cdot \vec{x} \, \, \varphi_j\,\varphi_i \; \mathrm{d}\bx
= 0.
\end{multline*}

On introduit alors les matrices suivantes, appelées matrices d'éléments finis. Les matrices $\bm M$ (de masse) et $\bm S$ (de rigidité) sont définies comme suit:
\begin{equation*}
\bm M_{kl} =\int_{\Omega_i} \varphi_k \varphi_l \, \mathrm{d}\bx, \quad \bm S_{kl} =\int_{\Omega_i} D_{i}\, \nabla \varphi_k \cdot \nabla \varphi_l \, \mathrm{d}\bx. \quad \begin{dcases}
k = 1,\dots,N_P, \\ l = 1,\dots,N_P,
\end{dcases}
\end{equation*}

La matrice $\bm J$ est similaire à la matrice de masse, mais elle est mise à l'échelle avec le coefficient $\vec{g}\cdot \vec{x}$:
\begin{equation*}
\bm J_{kl} =\int_{\Omega_i} \vec{g}\cdot \vec{x}\, \varphi_k \varphi_l \, \mathrm{d}\bx, \quad k = 1,\dots,N_P, \; l = 1,\dots,N_P,
\end{equation*}

Les coefficients dépendant du temps dans la base d'éléments finis de l'approximation de la magnétisation $M_i$ satisfont le système d'ODE suivant:
\begin{equation*}%\label{eq:BTPDE_Matrixform}
\bm M \frac{\partial U}{\partial t}=-\Bigl(\imunit\gamma f(t)\,{\bm J}+\bm S\Bigl)\, U\quad \text{où} \quad U = \left[\begin{array}{c}U_1\\ \vdots \\ U_{N_P}\end{array}\right].
\end{equation*}
Une méthode de résolution d'EDO du type Runge-Kutta est utilisée pour résoudre ce système d'équations semi-discrétisé.

Pour obtenir une approximation de la magnétisation sur l'ensemble du domaine $\Omega$, nous considérons tous les compartiments $\Omega_i$, avec $i\in \{1, \dots, N_\textup{cmpt}\}$, ensemble.

Chaque maillage d'éléments finis est alors constitué de
\begin{enumerate}
\item une liste de $N_\textup{element}$ éléments tétraédriques,
\item une liste de $N_\textup{node}$ n{\oe}uds:
\[
(\vec{Q}_1, \dots, \vec{Q}_{N_\textup{node}}) \in \mathbb{R}^{3 \times N_\textup{node}}.
\]
Un~n{\oe}ud correspond à un couple $(P,i)$ où $P$ est le sommet d'un tétraèdre~$T_l$ contenu dans $\Omega_i$. Ainsi un sommet $P$ apparaît $n_P$ fois dans la liste des n{\oe}uds, où $n_P$ est le nombre de $i\in\{1,\ldots, N_\textup{cmpt}\}$ tels que $P\in T_l\subset \Omega_i$.
\end{enumerate}

Les fonctions de base définies précédemment dans un compartiment~$\Omega_i$ fixé sont prolongées par $0$ en dehors de $\Omega_i$, et après renumérotation selon la liste des n{\oe}uds $Q_k$, on obtient une base de fonctions, appelée base nodale, sur tout $\Omega$:
\begin{equation*}
\varphi_k: \Omega \to [0, 1], \quad k \in \{1, \dots, N_\textup{node}\}.
\end{equation*}

Toute fonction $u$ dans l'espace des éléments finis peut être écrite comme une combinaison linéaire des fonctions de la base nodale ci-dessus:
\begin{equation*}
u(\vec{x},t) = \sum_{k = 1}^{N_\textup{node}} U_k(t) \varphi_k(\vec{x}) = \vec{U}^\transpose \vec{\varphi}(\vec{x}),
\end{equation*}
où $\vec{U} = (U_1(t), \dots, U_{N_\textup{node}}(t))^\transpose$ est le vecteur des coefficients et $\vec{\varphi} = (\varphi_1, \dots, \varphi_{N_\textup{node}})^\transpose$ est le vecteur des fonctions de base nodales des éléments finis.

\subsection{Représentation analytique de la solution par des fonctions propres de Laplace}

Dans cette partie, nous utilisons le formalisme matriciel de \cite{Callaghan1997,Barzykin1999} pour exprimer la magnétisation finale et le signal IRM de diffusion pour la séquence PGSE à l'aide des fonctions et valeurs propres de l'opérateur de Laplace généralisé\footnote{Ici, l'opérateur $\nabla \cdot D \nabla$ avec les conditions d'interface perméable décrites sera appelé opérateur de Laplace généralisé, par opposition à l'opérateur de Laplace pur $\nabla \cdot \nabla$.} sur le domaine connexe $\Omega = \bigcup_{i=1}^{N_\textup{cmpt}} \Omega_i$.

On dit que $\phi\in L^2(\Omega)$ est fonction propre normalisée de l'opérateur de Laplace généralisé associée à la valeur propre $\lambda$ si sa restriction~$\phi^i(\vec{x})$ au compartiment $\Omega_i$, pour $i\in\{1,\ldots, N_\textup{cmpt}\}$ satisfait~à
\begin{equation} \label{eq:laplace_eig}
-\nabla \cdot D_i \nabla \phi^i(\vec{x}) = \lambda \phi^i(\vec{x}), \quad \vec{x} \in \Omega_i, \ i \in \{1, \dots, N_\textup{cmpt}\},
\end{equation}
On demande de plus que $\phi^i$ satisfasse aux mêmes conditions de limi\-te et d'interface que pour l'équation de Bloch-Torrey, c'est-à-dire (on~rappelle qu'on a supposé que les membranes sont imperméables et donc, il n'y a pas d'interactions entre les domaines $\Omega_i$):
\begin{equation}\label{eq:lap_neumann}
D_i \nabla \phi^i(\vec{x}) \cdot \vec{n}_i(\vec{x}) = 0, \quad\vec{x} \in \partial \Omega_i, \ i \in \{1, \dots, N_\textup{cmpt}\}.
\end{equation}
On note $\left\{(\phi_n, \lambda_n)\right\}_{n \in \mathbb{N}^*}$ les solutions des équations ci-dessus \eqref{eq:laplace_eig}--\eqref{eq:lap_neumann}. Nous supposons que les valeurs propres sont réelles, positives ou nulles et ordonnées de manière croissante:
\begin{equation*}
0 = \lambda_1 \soutnew{<}{\leq} \lambda_2 \leq \lambda_3 \leq \cdots
\end{equation*}
Si le domaine $\Omega_i$ est connexe, seule la première valeur propre est nulle, et la fonction propre correspondante
est une fonction constante.

Pour ce qui nous intéresse, nous allons étudier la projection orthogonale de la magnétisation sur l'espace engendré par les $N_\textup{eig}$ premières fonctions propres $\phi_i$, où $N_\textup{eig}$ est un entier fixé en fonction de la géométrie de la cellule et la physique du système. Le choix de~$N_\textup{eig}$ est expliqué dans la section \ref{sect-echelle} ci-dessous.

Soit $\mat{L}$ la matrice diagonale, dont la diagonale est formée des $N_\textup{eig}$ premières valeurs propres de l'opérateur généralisé de Laplace:
\begin{equation*}%\label{eq:L_matrix}
\mat{L} = \operatorname{diag}(\lambda_1, \lambda_2, \dots,\lambda_{N_\textup{eig}}) \in \mathbb{R}^{N_\textup{eig} \times N_\textup{eig}}.
\end{equation*}

Soit $\mat{A}(\vec{g})$ la matrice $N_\textup{eig} \times N_\textup{eig}$ définie par:
\begin{equation*}
\mat{A}(\vec{g}) = g_x \mat{A}^x + g_y \mat{A}^y + g_z \mat{A}^z,
\end{equation*}
où $\vec{g} = (g_x, g_y, g_z)^\transpose$ est le vecteur gradient et $\mat{A}^x$, $\mat{A}^y$, et $\mat{A}^z$ sont trois matrices symétriques $N_\textup{eig} \times N_\textup{eig}$ dont les entrées sont les moments de premier ordre dans les directions des coordonnées du produit des paires de fonctions propres:
\begin{align*}
A_{mn}^x & = \int_{\Omega_i} x \, \phi_m(\vec{x}) \phi_n(\vec{x}) \, \mathrm{d}\bx, \quad & & (m, n) \in \{1, \dots, N_\textup{eig}\}^2, %\label{eq:lapeigfirstmm_x}
\\
A_{mn}^y & = \int_{\Omega_i} y \, \phi_m(\vec{x}) \phi_n(\vec{x}) \, \mathrm{d}\bx, \quad & & (m, n) \in \{1, \dots, N_\textup{eig}\}^2, %\label{eq:lapeigfirstmm_y}
\\
A_{mn}^z & = \int_{\Omega_i} z \, \phi_m(\vec{x}) \phi_n(\vec{x}) \, \mathrm{d}\bx, \quad & & (m, n) \in \{1, \dots, N_\textup{eig}\}^2. %\label{eq:lapeigfirstmm_z}
\end{align*}
Alors l'opérateur de Bloch-Torrey
\begin{equation*}
-\nabla \cdot D \nabla + \imunit \gamma \vec{g} \cdot \vec{x}
\end{equation*}
dans la base formée des $N_\textup{eig}$ premières fonctions propres de l'opérateur de Laplace généralisé est donné par la matrice $\mat{K}(\vec{g})$ à valeurs complexes:
\begin{equation*}
\mat{K}(\vec{g}) = \mat{L} + \imunit \gamma \mat{A}(\vec{g}).
\end{equation*}
La matrice suivante
\begin{equation*}%\label{eq:MF}
\mat{H}(\vec{g}, f) = \euler^{-\delta \mat{K}^*} \euler^{-(\Delta - \delta) \mat{L}} \euler^{-\delta \mat{K}},
\end{equation*}
donne la représentation du formalisme matriciel de la solution de l'équation de Bloch-Torrey pour la séquence PGSE.

En notant $\vec{\phi} = (\phi_1, \dots, \phi_{N_\textup{eig}})^\transpose$ le vecteur des fonctions propres de Laplace et $\vec{\nu} = \int_\Omega \rho(\vec{x}) \vec{\phi}(\vec{x}) \, \mathrm{d}\bx$ les coefficients de la magnétisation initiale dans la base de la fonction propre de Laplace, la projection de la magnétisation à la fin de la séquence de gradient est donnée par
\begin{equation*}
M^\text{MF}(\vec{x}, T_\textup{e}) = \vec{\phi}^\transpose\!(\vec{x}) \mat{H}(\vec{g},f) \vec{\nu}.
\end{equation*}

Le signal correspondant est
\begin{equation*}%\label{eq:mf_signal}
S^\text{MF}(\vec{g}, f) = \soutnew{\sqrt{|\Omega|} \mat{H}_1(\vec{g}, f) \vec{\nu}}{\vec{\Phi}^\transpose \mat{H}(\vec{g}, f) \vec{\nu}},
\end{equation*}
\soutnew{où $\mat{H}_1$ désigne la première ligne de $\mat{H}$ et $\int_\Omega \vec{\phi}(\vec{x}) \, \mathrm{d}\bx = (\sqrt{||\Omega|}, 0, \dots, 0)^\transpose$. Ceci est dû au fait que seule la première composante de Laplace reste après intégration, puisque la première fonction propre est constante, et que les autres fonctions propres sont orthogonales à $\phi_1$}{où $\vec{\Phi} = \int_\Omega \vec{\phi}(\vec{x}) \, \mathrm{d}\bx$. Si la densité de spin $\rho$ est uniforme dans tous les compartiments ($\rho_i = \rho$ pour tous les $i$), alors $\vec{\nu} = \rho \vec{\Phi}$. Si, en plus, tous les compartiments sont reliés au même domaine par des membranes perméables, l'expression se simplifie en $S^\text{MF}(\vec{g},f) = \rho |\Omega| H_{1 1}(\vec{g}, f)$, où $H_{1 1}$ est l'élément de la première ligne et de la première colonne de $\mat{H}$, car seule la fonction propre constante reste après intégration ($\vec{\Phi} = (\sqrt{|\Omega|}, 0, \dots, 0)^\transpose$)}.

Une quantité importante qui peut être dérivée du signal IRM de diffusion est
le \og coefficient de diffusion apparent\fg (\og apparent diffusion coefficient\fg), noté ADC et défini comme suit
(on suppose que le profil temporel $f$ est fixe, et on considère alors le signal comme une fonction de la valeur de $b$):
\begin{equation*}%\label{eq:ADCdef}
\ADC:= \left. -\frac{\partial}{\partial b} \log{\frac{S(b)}{S(0)}}\right\vert_{b=0}.
\end{equation*}
L'ADC donne une indication de la distance moyenne quadratique parcourue par les molécules d'eau dans la direction du gradient,
$\bug = \vec{g}/\|\vec{g}\|$,
moyennée sur toutes les positions de départ.

À partir du signal du formalisme matriciel, l'expression analytique de l'ADC est la suivante
\soutnew{}{(pour des applications, voir \cite{Haddar2018})}:
\begin{equation*}
%\label{eq:DeffEigenFunc}
\frac{\ADC}{\Dintr} = \sum_{n=1}^{N_\textup{eig}} \frac{(\bug \cdot \ba_{1n})^2\lambda_n\int_0^{T_E} F(t)
\left(\int_0^t e^{-\lambda_n(t - s)}f(s)\,\mathrm{d}s \right)\mathrm{d}t}{\Dintr \int_0^{T_E}F^2(t)\,\mathrm{d}t},
\end{equation*}
où \soutnew{}$F(t) \equiv \int_0^t f(s) ds$ et $\ba_{1n}=\bigl[ a_{1n}^x, a_{1n}^y, a_{1n}^z \bigr]^\transpose$ est la \og direction de diffusion\fg
du $n$-ième mode propre, dont les composantes sont les premiers moments
dans les 3 directions des axes principaux de la fonction propre associée, c'est-à-dire:
\begin{equation*}
\begin{split}
a_{1n}^x & = \frac{1}{\sqrt{\vert\Omega\vert}}\int_{\Omega} x \phi_n(\bx) \mathrm{d}\bx,\;\\
a_{1n}^y &= \frac{1}{\sqrt{\vert\Omega\vert}}\int_{\Omega} y \phi_n(\bx) \mathrm{d}\bx,\;\\
a_{1n}^z & = \frac{1}{\sqrt{\vert\Omega\vert}}\int_{\Omega} z \phi_n(\bx) \mathrm{d}\bx.
\end{split}
\end{equation*}
Nous rappelons au lecteur que la première fonction propre de Laplace est la fonction constante.

Pour clarifier la relation entre le ADC et la direction de codage de diffusion $\bug$,
nous réécrivons le ADC du formalisme matriciel comme suit:
\begin{equation*}
\frac{\ADC(\bug,f)}{\Dintr} = \bug^\transpose{}\,\frac{D^{\text{MF}}(f)}{\Dintr}\, \bug,
\end{equation*}
où le tenseur de diffusion effectif du formalisme matriciel est donné par:
\begin{equation*}
%\label{eq:dmf}
\frac{D^{\text{MF}}(f)}{\Dintr}= \sum_{n=1}^{N_\textup{eig}} J(\lambda_n,f)\, \ba_{1n}\, \ba_{1n}^\transpose,
\end{equation*}
avec $J(\lambda_n,f)$ dépendant de $\lambda_n$ et de $f$:
\begin{equation*}
J(\lambda_n,f)
=\frac{\lambda_n\int_0^{T_E} F(t)\left(\int_0^t e^{-\lambda_n(t-s)}f(s)\,\mathrm{d}s \right)\mathrm{d}t}{\Dintr\int_0^{T_E}F^2(t)\,\mathrm{d}t}.
\end{equation*}

Pour caractériser la contribution directionnelle du mode propre, nous utilisons le fait que
sa contribution dans la direction $\bug$ est $ J(\lambda_n,f)(\bug \cdot \ba_{1n})^2$.

\subsection{Échelle de longueur de la fonction propre}\label{sect-echelle}

La méthode des éléments finis développée dans la section \ref{sect-elts finis}, permet d'obtenir une approximation des valeurs propres $\lambda_n$ et des fonctions propres $\phi_n$, $1\leq n\leq N_\textup{node}$ (où $N_\textup{node}$ est le nombre de n{\oe}uds d'éléments finis) en fonction des éléments de la base nodale d'éléments finis $\varphi_k$, \hbox{$1\leq k\leq N_\textup{node}$}. Ceci est détaillé dans la section suivante. On note $\vec{p}_n$ les coordonnées de l'approximation de $\phi_n$ dans la base nodale d'éléments finis.

Notre objectif est de comparer les modes propres du signal calculé par le formalisme matriciel avec ceux du signal de référence issu de la résolution de l'EDP de Bloch-Torrey.

Nous ne voulons pas utiliser l'ensemble des valeurs propres et des vecteurs propres $\left\{(\lambda_n, \vec{p}_n)\right\}_{1 \leq n \leq N_\textup{node}}$. En effet, $N_\textup{node}$ est très grand (de l'ordre de $100\; 000$) et la plupart des modes propres à oscillation rapide sont liés à la discrétisation par éléments finis, et non à la physique du problème. Pour établir un lien avec la physique de la diffusion dans la géométrie de la cellule, nous fixons un intervalle restreint dans lequel nous conservons les valeurs propres calculées. Pour cela, nous allons définir une échelle de longueur associée à une valeur propre, et nous fixerons l'intervalle à $[0, (\pi / L_{\min})^2 D]$, où $L_{\min}$ est la plus courte échelle de longueur d'intérêt dans la géométrie de la cellule. De cette façon, le nombre de modes propres calculés, $N_\textup{eig}$, sera beaucoup plus petit que $N_\textup{node}$.

L'idée principale pour définir l'échelle de longueur d'une valeur propre (ou d'une fonction propre) est de considérer tout d'abord les fonctions propres de l'opérateur de Laplace généralisé en dimension~$1$ sur un segment de droite de longueur $L$, soumises à des conditions de Neumann nulles.

Sur un segment $[0,L]$ et pour une diffusivité $D$, de telles fonctions propres sont les fonctions $\psi$ définies sur $[0,L]$ telles que
\[
-D\psi''(x)=\mu \psi(x) \text{ et } \psi'(0)=\psi'(L)=0.
\]
Lorsqu'on impose $\mu\geq 0$, elles sont données, à constante multiplicative près, par
\[
\psi_n(x)=\cos \Bigl(\sqrt{\sfrac{\mu_n}{D}}\,x \Bigr)\quad\text{avec }\quad \mu_n=D\Bigl(\frac{(n-1)\pi}{L}\Bigr)^2,\; n\in\mathbb N^*.
\]
car $\psi_n'(x)=-\sqrt{\sfrac{\mu_n}{D}}\,\sin \bigl(\sqrt{\sfrac{\mu_n}{D}}\,x \bigr)$ et donc $\psi_n'(L)=0$ si et seulement si $\sqrt{\sfrac{\mu_n}{D}}\,L\in \pi\mathbb N$.

Pour $n>1$, la fonction $\psi_n$ a pour période $\spfrac{2L}{n-1}$ et donc pour fréquence $\psfrac{n-1}{2L}$.

En excluant la fréquence nulle, nous pouvons définir une coupure de fréquence maximale en choisissant $L_{\min}$ qui se traduira par une coupure pour les valeurs propres que nous calculons:
\begin{equation*}
0 \leq \mu_n \leq (\pi / L_{\min})^2 D.
\end{equation*}

En dimension $3$, nous appliquons la même coupure pour les valeurs propres:
\begin{equation*}
0 \leq \lambda_n \leq (\pi / L_{\min})^2 D.
\end{equation*}
Pour faire le lien entre la valeur propre calculée et l'échelle de longueur physique pertinente du problème, nous associons une échelle de longueur à une valeur propre en posant:
\begin{equation*}%\label{eq:lap_eig_lscale}
L(\lambda) =
\begin{dcases}
\soutnew{}{+ \infty, \quad} & \soutnew{}{\lambda = 0.} \\
\pi \sqrt{\sfrac{D}{\lambda}}, \quad & \lambda > 0,
\end{dcases}
\end{equation*}
et caractérisons la fonction propre calculée (appelée également mode propre) par $L(\lambda_n)$ au lieu de $\lambda_n$.

L'échelle de longueur associée à un mode propre est également appelée échelle spatiale du mode propre.

\subsection{Discrétisation de l'opérateur de Laplace par éléments finis}

Pour discrétiser l'opérateur de Laplace avec des conditions aux limites imperméables sur les membranes des cellules biologiques, nous construisons, comme dans la section \ref{sect-elts finis}, les matrices d'éléments finis suivantes: $\mat{M}, \mat{S}\in \mathbb{R}^{N_\textup{node} \times N_\textup{node}}$, connues dans la littérature comme les matrices de masse et de rigidité.

La discrétisation par éléments finis décrite ci-dessus transforme le problème continu des valeurs propres de l'opérateur de Laplace \eqref{eq:laplace_eig} en un problème discret et généralisé des valeurs propres de matrices: trouver $(\lambda, \vec{p}) \in \mathbb{R} \times \mathbb{R}^{N_\textup{node}}$ tel que\vspace*{-3pt}
\begin{equation} \label{eq:matrixeig}
\lambda \mat{M} \vec{p} = \mat{S} \vec{p},
\end{equation}
dont on retiendra les $N_\textup{eig}$ plus petites valeurs propres et les vecteurs propres correspondants $\{(\lambda_n, \vec{p}_n)\}_{1 \leq n \leq N_\textup{eig}}$, avec $N_\textup{eig} \leq N_\textup{node}$. Notons cependant qu'il existe au total $N_\textup{node}$ solutions au problème \eqref{eq:matrixeig}. En revenant à l'espace des éléments finis, l'approximation de la fonction propre $\phi_n(\vec{x})$, que l'on notera par la même lettre, associée à la valeur propre $\lambda_n$ s'écrit alors\vspace*{-3pt}
\begin{equation*}
\phi_n(\vec{x}) = \sum_{k = 1}^{N_\textup{node}} p_n^k \varphi_k(\vec{x}) = \vec{p}_n^\transpose \vec{\varphi}(\vec{x}), \quad n \in \{1, \dots, N_\textup{eig}\},
\end{equation*}
où les entrées du vecteur propre $\vec{p}_n$ sont les coefficients de la fonction propre $\phi_n$ dans la base des éléments finis.

En utilisant la notation matricielle, cette conversion peut également s'écrire $\vec{\phi} = \mat{P}^\transpose \vec{\varphi}$, où\vspace*{-3pt}
\begin{gather*}
\vec{\varphi} = (\varphi_1, \dots, \varphi_{N_\textup{node}})^\transpose,\ \vec{\phi} = (\phi_1, \dots, \phi_{N_\textup{eig}})^\transpose, \\
\tag*{et} \mat{P} = (\vec{p}_1, \dots, \vec{p}_{N_\textup{eig}}) \in \mathbb{R}^{N_\textup{node} \times N_\textup{eig}}.
\end{gather*}
Les intégrales des fonctions propres discrétisées par éléments finis sont alors données par\vspace*{-3pt}
\[
\vec{\Phi} = \int_\Omega \vec{\phi}(\vec{x}) \, \mathrm{d}\bx = \mat{P}^\transpose \mat{M} \vec{o},\quad \text{où}\quad\vec{o} = (1, \dots, 1)^\transpose \in \mathbb{R}^{N_\textup{node}}.
\]
De même, les coefficients de la densité de spin initiale dans la base de la fonction propre discrétisée par éléments finis sont donnés par\vspace*{-3pt}
\[
\vec{\nu} = \int_\Omega \rho(\vec{x}) \vec{\phi}(\vec{x}) \, \mathrm{d}\bx = \mat{P}^\transpose \mat{M} \vec{\rho},
\]
où $\vec{\rho} = (\rho_{i(k)})_{1 \leq k \leq N_\textup{node}} \in \mathbb{R}^{N_\textup{node}}$ et $i(k) \in \{1, \dots, N_\textup{cmpt}\}$ est tel que $k \in \mathcal{I}_{i(k)}$.

\section{Simulations numériques}

\label{Résultats}

Les résultats numériques de cette section concernent le neurone pyra\-midal \texttt{02b\_pyramidal1aACC},
dont le périmètre est
\[
[-70,113]\lunit \times[-197,165]\lunit\times[-14,18]\lunit.
\]
Le maillage par éléments finis de ce neurone comporte 44908 nœuds et 171017 éléments tétraédriques.

\subsection{Comparaison du signal de formalisme matriciel avec le signal de référence}

Dans les simulations ci-dessous, nous avons fixé la diffusivité intrinsèque en posant $\Dintr = 2\times 10^{-3} \dunit$.
Nous avons simulé deux séquences de codage de diffusion:
\begin{align*}
&\text{SEQ1 (PGSE, $\delta=10.6\tunit, \Delta=13\tunit$);}\\
&\text{SEQ2 (PGSE, $\delta=10.6\tunit, \Delta=73\tunit$).}
\end{align*}
L'ensemble des valeurs de $b$ simulées sont
\[
\{0,1000,2000,3000,4000\}\bunit.
\]
La géométrie est celle du neurone pyramidal.

Nous fixons l'échelle de longueur d'intérêt minimale à $l_s^{\min} = 4 \lunit$, de sorte que nous avons calculé les valeurs propres appartenant
à $[0, (\pi/l_s^{\min})^2 \Dintr]$.
Dans cet intervalle, $336$ valeurs propres (correspondant donc à $336$ fonctions propres de Laplace) ont été trouvées, dont $\lambda_1=0$, qui correspond à l'échelle de longueur $l_s = \infty$.
Il existe 6 modes propres avec une échelle de longueur $l_s > 100 \lunit$, ils correspondent aux échelles de longueur (arrondies au $\lunit$),
\[
\{405, 343, 162, 156, 133, 127, 106\}\lunit,
\]
respectivement.

Ensuite, nous comparons les tracés des signaux d'IRM de diffusion calculés de deux manières différentes:
\begin{enumerate}
\item Signaux de référence issus de la résolution de l'EDP de Bloch-Torrey \soutnew{}{($S^{\text{BTPDE}}$)}, calculés dans 151 directions de diffusion uniformément réparties dans la sphère unité;
\item Signaux issus du formalisme matriciel \soutnew{}{($S^{\text{MF}}$)} en utilisant les 336 fonctions propres trouvées dans l'intervalle $[0, (\pi/l_s^{\min})^2 \Dintr]$, $l_s^{\min}=4\lunit$, calculées dans 151 directions de diffusion uniformément distribuées dans la sphère unité;
\end{enumerate}

Nous voyons sur la figure \ref{fig:sighardi}, qu'à $b = 1000\bunit$,
les formes des signaux sont ellipsoïdes.
À $b = 4000\bunit$, les formes des signaux ne sont plus ellipsoïdes.

\begin{figure}[htb]
\centering
\includegraphics[width=0.22\textwidth]{btpde_151_ex1_b1000}
\includegraphics[width=0.22\textwidth]{mf_151_ex1_b1000}
\includegraphics[width=0.22\textwidth]{btpde_151_ex1_b4000}
\includegraphics[width=0.22\textwidth]{mf_151_ex1_b4000} \\
\includegraphics[width=0.22\textwidth]{btpde_151_ex2_b1000}
\includegraphics[width=0.22\textwidth]{mf_151_ex2_b1000}
\includegraphics[width=0.22\textwidth]{btpde_151_ex2_b4000}
\includegraphics[width=0.22\textwidth]{mf_151_ex2_b4000}
\caption{\label{fig:sighardi} Signaux BTPDE, $S^{\text{BTPDE}}/S_0$, en 151
directions de codage de diffusion. Signaux MF, $S^{\text{MF}}/S_0$, dans 151 directions de codage de diffusion.
Les points noirs indiquent la direction du codage par diffusion, multipliée par l'ampleur de l'atténuation du signal. La couleur indique la valeur de l'atténuation du signal. SEQ1 est (PGSE, $\delta=10.6\tunit, \Delta=13\tunit$), SEQ2 est (PGSE, $\delta=10.6\tunit, \Delta=73\tunit$).
La géométrie est celle du neurone pyramidal \texttt{02b\_pyramidal1aACC}. }
\end{figure}

\Subsection{La contribution de chaque mode propre au signal}

Comme $S^{\text{MF}}$ contient les contributions de tous les modes propres calculés dans l'intervalle demandé, pour avoir une idée de
l'importance de chaque mode propre, nous avons calculé la différence de signal qui résulte de la suppression d'un mode propre,
par rapport à l'utilisation de l'ensemble complet des modes propres calculés. Cette différence de signal est calculée pour chaque
valeur de $b$, en moyenne sur 30 directions de gradient, ces directions étant uniformément distribuées sur la sphère unité.
Pour la fonction propre $i$, la différence de signal est obtenue comme suit:
\begin{equation*}
E^{\text{RM}, i}(f,b) = \frac{\sum_{j=1}^{30}\left( S^\text{MF}(f,\bg_j)-S^{\text{MF}, \text{RM}, i }(f,\bg_j)\right)^2}{\sum_{j=1}^{30}\left( S^\text{MF}(f,\bg_j)\right)^2}.
\end{equation*}
Le signal $S^\text{MF}$ utilise l'ensemble complet des fonctions propres calculées, le signal $S^{\text{MF}, \text{RM}, i }$ exclut la $i$-ième fonction propre. Dans ce qui suit, les différences de signal seront données
pour deux séquences et à 2 valeurs de $b$, dans l'ordre suivant:
\begin{multline*}
\text{\{(SEQ1, $b = 1000\bunit$), (SEQ1, $b=4000\bunit$),}\\
\text{(SEQ2, $b=1000\bunit$), (SEQ2, $b=4000\bunit$)\}.}
\end{multline*}
Nous nous attendons à ce que la deuxième valeur soit la plus élevée et la troisième la plus basse.
Nous dirons que la $i$-ième fonction propre est \og significative\fg si $E^{\text{RM}, i}(f,b)$ est supérieur à un certain seuil.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.48\textwidth]{sigdiff_experi1_0-50} \hfill
\includegraphics[width=0.48\textwidth]{sigdiff_experi1_50-350}\\[5pt]
(a)\enspace SEQ1 (PGSE, $\delta=10.6\tunit, \Delta=13\tunit$);

\vspace{0.5 cm}

\includegraphics[width=0.48\textwidth]{sigdiff_experi3_0-50} \hfill
\includegraphics[width=0.48\textwidth]{sigdiff_experi3_50-350}\\[5pt]
(b)\enspace SEQ2 (PGSE, $\delta=10.6\tunit, \Delta=73\tunit$);
\end{center}
\caption{\label{fig:sigdiff}
Les différences de signal dues à la suppression de chaque mode propre, par rapport à
l'utilisation de l'ensemble complet des 336 modes propres.
Les valeurs propres ont été converties en une échelle de longueur.
La couleur indique la \og direction de diffusion\fg des modes propres, en fonction des valeurs du vecteur RVB $\bbd{c}_n$ qui est lié à $\ba_{1n}$.
La géométrie est celle du neurone pyramidal\\\texttt{02b\_pyramidal1aACC}.}
\end{figure}

Dans la figure \ref{fig:sigdiff}, nous montrons les modes propres significatifs.
Pour visualiser la \og direction de diffusion\fg des modes propres,
nous utilisons une échelle de couleurs RVB (rouge, vert, bleu) basée sur les valeurs de l'indice de gris. Ceci définit un
vecteur, appelé vecteur RVB et noté $\bbd{c}_n$ avec trois composantes positives ou nulles:\vspace*{-3pt}\enlargethispage{\baselineskip}%
\begin{equation*}
\bbd{c}_{n}=\frac{\bigl[
|a_{1n}^x|, |a_{1n}^y|, |a_{1n}^z|
\bigr]}{\sqrt{(a_{1n}^x)^2+ (a_{1n}^y)^2+(a_{1n}^z)^2}}.
\end{equation*}
(On rappelle que le vecteur $\ba_{1n}=\bigl[ a_{1n}^x, a_{1n}^y, a_{1n}^z \bigr]^\transpose$ est la \og direction de diffusion\fg du $n$-ième mode propre).

La couleur indiquée par le vecteur RVB $\bbd{c}_n$ peut être utilisée pour évaluer la contribution relative du mode propre dans les 3 principales directions de diffusion, $x$, $y$ ou $z$.

Pour SEQ1, les modes propres significatifs pour $25\lunit \leq l_s \leq 50\lunit$ sont principalement verts, ce qui signifie qu'ils contribuent à la diffusion dans la direction $y$. Lorsque $0 \leq l_s \leq 25\lunit$, il y a beaucoup plus de modes propres significatifs
qui sont rouges, ce qui signifie qu'ils contribuent principalement à la diffusion dans la direction $x$.
Il~n'y a qu'un seul mode qui est bleu, ce qui signifie qu'il contribue de manière significative à la diffusion dans la
direction $z$. Ces résultats étaient prévisibles puisque ce neurone se situe principalement dans le plan $x-y$.

Nous voyons également qu'avec une valeur de $b$ plus élevée, il y a plus de modes propres significatifs.
et les différences de signaux sont également plus élevées,
par rapport à une valeur de $b$ inférieure.
Le mode propre le plus significatif est celui dont l'échelle de longueur est $l_s(\lambda_n) = 33\lunit$, il est le plus aligné sur la direction $y$.
La~figu\-re~\ref{fig:lapeigs} montre cette fonction propre et nous voyons que l'échelle de longueur correspond à la \og longueur d'onde\fg des oscillations signi\-fi\-ca\-tives de la fonction propre dans la géométrie de la cellule.
Parmi les modes propres avec des échelles de longueur plus grandes, $ l_s > 50 \lunit$, il y a beaucoup moins de
modes propres significatifs qu'entre $0 \leq l_s \leq 50\lunit $
et ils sont principalement dans la direction $y$ (étant principalement verts).
Le mode propre correspondant à $l_s(\lambda_n) = 343.6\lunit$ est représenté sur la figure \ref{fig:lapeigs}.
Sa suppression entraînera des différences de signal de
\{19.1\%, 83.8\%, 8.8\%, 40.8\%\}. On peut voir que la \og longueur d'onde\fg de ce mode est
plus grande (oscillation plus lente) que celle du mode avec $l_s(\lambda_n) = 33\lunit$.

Pour des échelles de longueur petites, proches de $l_s = 15\lunit$, il~appa\-raît que beaucoup des modes propres significatifs pour SEQ1 ne le sont plus pour SEQ2.

Dans la figure \ref{fig:lapeigs}, nous montrons que le mode propre
bleu lorsque $l_s = 15.6\lunit$ est significatif à la fois pour SEQ1 et SEQ2.
Ce mode propre est principalement dans la direction $z$, et les oscillations rapides se trouvent dans les branches de la dendrite.
La fonction propre qui possède la plus longue échelle de $l_s = 405\lunit$, représentée sur la figure \ref{fig:lapeigs}, entraîne les différences de signal suivantes lorsqu'il est supprimé: $E^{\text{RM}, i}=\{1.06\%, 6.58\%, 0.57\%, 1.22\%\}$.\enlargethispage{\baselineskip}

\begin{figure}[htb]
\vspace*{-.5\baselineskip}
\noindent
\begin{tabular}{cc}
\textbf{(a)}&\textbf{(b)}\\
\includegraphics[width=0.45\textwidth]{Eigfun1}&
\includegraphics[width=0.45\textwidth]{Eigfun6}\\[15pt]
\textbf{(c)}&\textbf{(d)}\\
\includegraphics[width=0.45\textwidth]{Eigfun5}&
\includegraphics[width=0.45\textwidth]{eigfun_long1}
\end{tabular}
\caption{\label{fig:lapeigs} Les fonctions propres correspondant à l'échelle de longueur $l_s(\lambda_n) = 33\lunit$ \textbf{(a)}, $343.6\lunit$ \textbf{(b)},
$15.6\lunit$ \textbf{(c)} et $405\lunit$ \textbf{(d)}.
La géométrie est celle du neurone pyramidal.}
\vspace*{-.5\baselineskip}
\end{figure}


\section{Application à l'imagerie}

Pour étudier plus en profondeur le signal de diffusion IRM des neurones,
nous avons divisé les
neurones en composants géométriques disjoints: à savoir,
le soma et les branches de la dendrite.
À titre d'illustration, nous présentons sur la figure \ref{fig:neuron_broken} le neurone \texttt{03a\_spindle2aFI} en maillage tétraédrique volumique
divisé en sous-mailles du soma et des deux branches de la dendrite.
\begin{figure}[htb]
	\centering
	\begin{tikzpicture}
\draw (-2, 0) node[inner sep=0] {\includegraphics[height=6cm]{03a_spindle2aFI_soma_mesh.jpg}};
\draw (0.3, 0.5) node {\mbox{\larger[4]$+$}};
\draw (2.5, 0) node[inner sep=0] {\includegraphics[height=6cm]{03a_spindle2aFI_dendrites_mesh.jpg}};
\draw (5.0, 0.5) node {\mbox{\larger[4]$=$}};
\draw (7, 0) node[inner sep=0] {\includegraphics[height=6cm]{03a_spindle2aFI_FEmesh.jpg}};
	\end{tikzpicture}
	\caption{\label{fig:neuron_broken} Le maillage tétraédrique volumique du neurone \texttt{03a\_spindle2aFI} est
divisé en trois composantes géométriques déconnectées: le soma et deux branches de dendrites}
	\label{fig:03a_spindle2aFI_mesh}
\end{figure}

\subsection{Comportement à haute valeur de \texorpdfstring{$b$}{b}}
\label{sec:simul_highb}

Dans \cite{Veraart2019}, il a été montré expérimentalement que le signal IRM de diffusion des
structures tubulaires, telles que les axones, présentent un certain comportement à haute valeur de $b$.
Plus précisément, le signal moyen dans la direction de la diffusion, noté $S_\textup{ave}(b)$, est affine en $\sfrac{1}{\sqrt{b}}$ pour les valeurs élevées de $b$:
\begin{equation*}
%\label{eq:power_law}
S_\textup{ave}(b) \equiv \int_{\|\bug\| = 1} S_{\bug} (b) \,\mathrm{d}\bug \sim c_0 + c_1 \frac{1}{\sqrt{b}}.
\end{equation*}
Étant donné que les dendrites des neurones ont également une structure tubulaire, nous testons si
le signal moyen, $S_\textup{ave}(b)$, dans la direction de la diffusion des branches de dendrites montre également
le comportement ci-dessus à haute valeur de $b$. Nous avons calculé $S_\textup{ave}(b)$ pour le neurone entier ainsi que pour ses deux branches dendritiques,
correspondant à la moyenne sur 120 directions de gradient uniformément réparties
dans la sphère unité. Les résultats sont présentés dans la figure \ref{fig:high_bvalue_fitting}.
\begin{figure}[htb]
	\centering
	\includegraphics[width=0.49\textwidth]{fitting_sqrtb_4aACC_Dintri1e-03.png}
	\includegraphics[width=0.49\textwidth]{fitting_sqrtb_4aACC_Dintri2e-03.png}
	\caption{\Changelback
	Le signal moyen en direction de la diffusion pour le neurone \texttt{03b\_spindle4aACC}.
La moyenne de $S_\textup{ave}(b)$ est calculée sur 120 directions de diffusion, uniformément distribuées dans la sphère unité,
et il est normalisé de sorte que $S_\textup{ave}(b=0)=1$. Les paramètres de simulation sont\\
\mbox{}\hfill$rtol = 10^{-3}$, $atol = 10^{-5}$, $Htetgen = 0,5\mu m$.\hfill\mbox{}\\
La séquence de codage de la diffusion est\\
\mbox{}\hfill
PGSE ($\delta = 10\tunit$, $\Delta = 43\tunit$).\hfill\mbox{}\\
Les valeurs de $b$ sont\\
$b=\{60000, 40000, 20000, 12000, 10000, 8000, 7000$,\\
\mbox{}\hfill
$6000, 4000, 2500\} \bunit$.\\
\mbox{}\hfill
(a) $\Dintr = 1 \times 10^{-3}\dunit$. (b) $\Dintr = 2 \times 10^{-3}\dunit$.\hfill\mbox{}
}
	\label{fig:high_bvalue_fitting}
\end{figure}
On voit clairement la
relation affine entre $S_\textup{ave}(b)$ et $\sfrac{1}{\sqrt{b}}$
dans les branches de la dendrite pour des valeurs de~$b$ dans la gamme $2500\bunit \leq b \leq 20000\bunit$.
En simulant à la fois pour $\Dintr = 2\times 10^{-3} \dunit$ et $\Dintr = 1\times 10^{-3} \dunit$,
nous voyons que la pente ajustée $c_1$ est effectivement $\sfrac{1}{\sqrt{\Dintr}}$.

En revanche, dans le neurone entier, en raison de la présence du soma, cette relation affine n'existe pas.

\subsection{Biomarqueurs de la taille des somas}
\label{sec:biomarqueurs_soma}

Comme nous l'avons montré sur la figure \ref{fig:high_bvalue_fitting}, la relation affine entre $S_\textup{ave}(b)$ et $\sfrac{1}{\sqrt{b}}$, en~d'autres termes, la loi établie dans \cite{Veraart2019}, appelée loi de puissance dans la suite, ne tient pas en raison de la présence du soma et des effets d'échange entre le soma et les dendrites. En exploitant la collection de maillages de neurones réalistes, nous montrons statistiquement dans cette section que la déviation de la loi de puissance a le potentiel de servir de biomarqueur pour révéler la taille du soma.

\begin{figure}[htb]
	\centering
	\includegraphics[width=0.7\textwidth]{03a_spindle2aFI_inflection.png}
	\caption{
La courbe de signal, moyenné dans la direction, pour le neurone \texttt{03a\_spindle2aFI}.
Les signaux sont calculés numériquement à l'aide du module de formalisme matriciel de la boîte à outils SpinDoctor.
La moyenne de $S_\textup{ave}(b)$ a été calculée sur 64 directions de diffusion, uniformément distribuées dans la sphère unité, et elle est normalisée de sorte que $S_\textup{ave}(b=0)=1$. Les valeurs de $b$ sont supérieures à 278 $\bunit$ et la diffusivité est de $\Dintr = 2 \times 10^{-3}\dunit$.
L'amplitude du gradient est constante, $\gamma |\bg| = 10^{-5} \text{s}^{-1}\cdot \text{mm}^{-1}$, et on a fait varier $\delta$ pour obtenir une large gamme de valeurs de~$b$, tout en choisissant $\Delta=\delta$ (séquence PGSE).
	Le point bleu indique le point d'inflexion de la courbe du signal simulé. La loi de puissance est ajustée autour du point d'inflexion. La région de la loi de puissance est la largeur de la plage où l'erreur relative entre le signal simulé et l'approximation de la loi de puissance est inférieure à 2\%. La zone située entre la courbe simulée et la loi de puissance à gauche du point d'inflexion représente l'erreur d'approximation de la loi de puissance.}
	\label{fig:possible_biomarkers}
\end{figure}

Pour ce faire, nous avons réalisé les simulations suivantes qui sont légèrement différentes des expériences constantes ($\delta$,$\Delta$) établies dans
\cite{Veraart2019} et illustrées sur la figure \ref{fig:high_bvalue_fitting}. Dans ce qui suit, nous avons retenu l'amplitude du gradient
constante, $\gamma|\bg| = 10^{-5} s^{-1}\cdot mm^{-1}$, et fait varier $\delta$ pour obtenir une large gamme de valeurs de $b$, tout en choisissant $\Delta=\delta$ (séquence PGSE). Les simulations ont été effectuées dans 64 directions de gradient
et les signaux ont été moyennés sur ces directions. Cette opération a été réalisée pour un ensemble de 65 neurones.

Dans la figure \ref{fig:possible_biomarkers}, nous montrons un exemple de la courbe du signal simulé et de l'approximation de la loi de puissance pour le neurone \texttt{03a\_spindle2aFI}. À partir des signaux simulés moyennés en direction, nous trouvons le point d'inflexion (point bleu) de la courbe du signal (courbe noire). Nous ajustons la loi de puissance (ligne droite en pointillés bleus) autour du point d'inflexion.
\begin{figure}[htb]
	\centering
	\includegraphics[width=0.98\textwidth]{soma_volume_vs_biomarkers.png}
	\caption{(a) le logarithme du volume somatique en fonction de la coordonnée x du point d'inflexion $x_0$; (b) le logarithme du volume du soma en fonction de la coordonnée y du point d'inflexion $y_0$; (c) le logarithme du volume somatique en fonction de l'ordonnée à l'origine de la loi de puissance $c_0$; (d) le logarithme du volume du soma en fonction de la pente de la loi de puissance $c_1$; (e) le logarithme du volume somatique en fonction de l'erreur d'approximation de la loi de puissance $\mathcal{E}$; (f) le logarithme du volume du soma en fonction de la largeur de la région de la loi de puissance $w$. Chaque point bleu représente les données d'un des 49 neurones (28 neurones fusiformes et 21 neurones pyramidaux) retenus pour cette étude.}
	\label{fig:volume_biomarqueurs}
\end{figure}

La région de la loi de puissance est la plage où l'erreur relative entre la courbe du signal simulé et l'ajustement de la loi de puissance est inférieure à 2\% (largeur de la région jaune) et l'erreur \hbox{d'approximation} est estimée par la zone entre la courbe du signal et l'ajustement de la loi de puissance \emph{à gauche du point d'inflexion} (la zone verte).

Afin de caractériser l'influence du soma sur l'approximation de la loi de puissance, nous avons choisi les 6 biomarqueurs candidats suivants:\vspace*{-3pt}\enlargethispage{1.5\baselineskip}
\begin{itemize}
	\item $x_0$: la coordonnée x du point d'inflexion;
	\item $y_0$: la coordonnée y du point d'inflexion;
	\item $c_0$: l'ordonnée à l'origine de l'ajustement de la loi de puissance;
	\item $c_1$: la pente de l'ajustement de la loi de puissance;
	\item $\mathcal{E}$: l'erreur d'approximation de la loi de puissance;
	\item $w$: la largeur de la région de la loi de puissance.
\end{itemize}

Une étude statistique des 6 biomarqueurs candidats ci-dessus sur l'ensemble des 65 neurones a été réalisée. Comme le sous-échantillonnage lorsque $\sfrac{1}{\sqrt{b}}$ s'approche de $0$ pourrait produire une erreur numérique significative, nous n'avons gardé que les neurones dont les $x_0$ sont supérieurs à 0,016 $\text{mm}\cdot \text{s}^{-1/2}$. Au total, 28 neurones fusiformes et 21 neurones pyramidaux ont été retenus.

Nous traçons d'abord les biomarqueurs candidats par rapport au volume du soma $v_\textup{soma}$ dans la figure \ref{fig:volume_biomarqueurs}. Chaque point de données dans la figure correspond à un neurone (pour un total de 49). On constate que $x_0$, $c_0$, $c_1$, $\mathcal{E}$ et $w$ présentent une relation exponentielle avec le volume du soma. Les équations ajustées nous permettent de déduire le volume du soma en mesurant les biomarqueurs. Nous constatons également
que $y_0$ n'est pas un biomarqueur pour le volume du soma.

\begin{figure}[htb]
	\centering
	\includegraphics[width=0.98\textwidth]{soma_fraction_vs_biomarkers.png}
	\caption{(a) la fraction de volume du soma en fonction de la coordonnée x du point d'inflexion $x_0$; (b) fraction volumique du soma en fonction de la coordonnée y du point d'inflexion $y_0$; (c) la fraction volumique du soma en fonction de l'ordonnée à l'origine de la loi de puissance $c_0$; (d)~la fraction volumique du soma en fonction de la pente de la loi de puissance $c_1$; (e) la fraction volumique du soma en fonction de l'erreur d'approximation de la loi de puissance~$\mathcal{E}$; (f) fraction volumique du soma en fonction de la largeur de la région de la loi de puissance $w$. Chaque point bleu représente les données d'un des 49 neurones (28 neurones fusiformes et 21 neurones pyramidaux) retenus pour cette étude.}
	\label{fig:fraction_biomarqueurs}
\end{figure}

De même, nous montrons le diagramme de dispersion des biomarqueurs candidats par rapport à la fraction $f_\textup{soma}$ du volume du soma par rapport au volume total du neurone dans la figure \ref{fig:fraction_biomarqueurs}. Dans ce cas, les $x_0$, $c_1$ et $w$ ne sont pas des biomarqueurs de la fraction volumique du soma. Les biomarqueurs candidats $y_0$, $c_0$ et $\mathcal{E}$ semblent capables d'indiquer la limite inférieure de la fraction volumique du soma.\enlargethispage{-2\baselineskip}

Nous notons que l'objectif de cette section est de donner un exemple de la recherche possible qui peut être menée. Une étude plus systématique est nécessaire pour obtenir des biomarqueurs plausibles pour la taille du soma, mais cela n'entre pas dans le cadre de cet article.

\Changelback
\backmatter
\bibliographystyle{jepplain+eid}
\bibliography{xups22-03}
\end{document}

