# Modélisation macroscopique de la croissance des tumeurs cérébrales

Olivier Clatz, Pierre-Yves Bondiau, Hervé Delingette, Nicholas Ayache

► To cite this version:

Olivier Clatz, Pierre-Yves Bondiau, Hervé Delingette, Nicholas Ayache. Modélisation macroscopique de la croissance des tumeurs cérébrales. 17e Congrès Français de Mécanique, Aug 2005, Troyes, France. hal-04262501

HAL Id: hal-04262501

https://inria.hal.science/hal-04262501v1

Submitted on 27 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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 CC BY 4.0 - Attribution - International License
# Modélisation macroscopique de la croissance des tumeurs cérébrales.

O. Clatz*, P.Y. Bondiau, H. Delingette and N. Ayache

Projet de recherche Epidaure, INRIA Sophia Antipolis
2004 route des Lucioles, 06902 Sophia Antipolis, France

*Contact : Olivier.Clatz@sophia.inria.fr

## Résumé :

Nous proposons dans cet article une méthode permettant de simuler la croissance du glioblastome (GBM), la tumeur gliale (cérébrale) la plus agressive. Ce type de simulation trouve de nombreuses applications, à la fois en radiothérapie pour une meilleure quantification de la dose administrée, et en neurochirurgie dans l'aide à la planification de l'opération et la définition des marges de résection de la tumeur. La simulation repose sur une modélisation couplée de deux effets : l'invasion par diffusion du glioblastome et la déformation mécanique des structures cérébrales voisines. Le premier phénomène est modélisé par une équation de type réaction-diffusion alors que le second repose sur les lois de comportement de la mécanique des milieux continus. Une nouvelle loi de comportement est proposée, permettant de coupler localement les efforts mécaniques à la densité de cellules tumorales dans le parenchyme cérébral. Le problème est résolu avec la méthode des éléments finis. Cette résolution repose sur un maillage tétraédrique du cerveau spécifique au patient, construit de manière automatique par des outils d'analyse d'images médicales. Enfin, une première évaluation du modèle est effectuée en simulant la croissance d'un glioblastome chez un patient donné. Cette croissance in-silico est comparée à la croissance réellement observée dans l'IRM du même patient 6 mois plus tard.

## Abstract :

We propose a new model to simulate the growth of glioblastomas multiforma (GBM), the most aggressive glial tumors. Different medical applications could benefit from this kind of simulation, including an optimized dosimetry in radiotherapy or a better neuro-surgical planning in case of tumor resection. Because the GBM shows a preferential growth in the white fibers and have a distinct invasion speed with respect to the nature of the invaded tissue, we rely on an anatomical atlas to introduce this information into the model. This atlas includes a white fibers diffusion tensor information and the delineation of cerebral structures having a distinct response to the tumor aggression. We use the finite element method (FEM) to simulate both the invasion of the GBM in the brain parenchyma and its mechanical interaction (mass effect) with the invaded structures. In addition, we propose a new coupling equation taking into account the mechanical influence of the tumor cells on the invaded tissues. This tumor growth model is assessed by comparing the in-silico GBM growth with the real GBM growth observed between two magnetic resonance images (MRIs) of a patient acquired with six months difference.

# 1 Introduction

La majorité des tumeurs du système nerveux central sont d'origine gliale. Parmi celles-ci, les glioblastomes (GBMs) sont considérés comme les plus agressifs. Peu de travaux se sont intéressés au développement de modèles numériques macroscopiques spécifiques au patient permettant de prédire le développement de ces tumeurs. De tels modèles trouvent de nombreuses applications dans le traitement de la pathologie. En radio-thérapie, une bonne connaissance de l'envahissement spatial de la tumeur permettrait de prédire la densité de cellules tumorales dans le parenchyme cérébral, et plus spécifiquement dans les zones ou la sensibilité du scanner ne permet pas de détecter ces cellules. Les marges dosimétriques de radiothérapie, pour l'instant choisies à distance constante de la limite "visible" de la tumeur, pourraient ainsi être adaptées à la prédiction réalisée par le modèle. En neuro-chirurgie, la prédiction de l'envahissement de la tumeur dans les différentes aires cérébrales peut être corrélée avec un atlas fonctionnel permettant d'anticiper les déficits cognitifs liés à la pathologie. Cette donnée se révèle fondamentale dans le choix des zones de résection définies lors du planning pré-opératoire.

Dans cet article, nous proposons une évolution d'un modèle de propagation existant ainsi qu'une nouvelle loi de couplage entre la densité de cellules tumorales et les équations de la mécanique. A partir de l'image par résonance magnétique (IRM) du patient, nous utilisons différents outils classiques d'analyse d'images médicales
pour segmenter le cerveau du patient, ainsi que les sous structures ayant un comportement distinct vis à vis de la tumeur. Ces segmentations sont utilisées pour la construction et l'initialisation d'une modèle tetraédrique spécifique au patient. Pour la première fois les résultats issus de la simulations sont confrontés aux données enregistrées chez le même patient plusieurs mois plus tard.

Les travaux présentés dans cet article démontrent la faisabilité de la prédiction de l'envahissement du glioblastome chez un patient spécifique, avec une comparaison sur des données réelles. Ils représentent la première étape en vue de futures études à une plus grande échelle permettant une validation. Cependant, plusieurs valeurs sont ajustées manuellement et nécessiteraient des mesures in-vivo non disponibles pour leur identification. De plus, certaines de ces valeurs sont des paramètres intrinsèques de la tumeur et varient certainement d'une tumeur à l'autre.

# 2 Modèle macroscopique d'envahissement de la tumeur

## 2.1 Évolution de la densité de cellules tumorales

On peut classifier les principaux travaux relatifs à la modélisation de la croissance des tumeurs suivant l'échelle de modélisation. Au niveau microscopique, ces modèles décrivent le comportement individuel de chaque cellule tumorale, les interactions inter-cellules et avec leur environnement. Ce niveau de détail est souvent associé à une grande complexité numérique lorsque l'on s'intéresse à des phénomènes de croissance prolongés dans le temps (plusieurs mois) et impliquant un très grand nombre de cellules (automates cellulaires). À un niveau macroscopique, les modèles décrivent l'évolution de quantités globales, comme la densité de cellules tumorales. La majeure partie des travaux dans ce domaine utilisent une équation de type réaction-diffusion isotrope pour décrire l'évolution de cette densité [1, 2] :

$$ \frac{\partial c}{\partial t} = -div(J) + S(c, t) \quad (1) $$

- c représente la densité de cellules tumorales.

- J représente le flux de cellules tumorales par diffusion<fnref n="1" />, obéissant à la loi de Fick : $J = -d \nabla c$, avec d coefficient de diffusion.

- $S(c, t)$ représente le facteur source. Classiquement $S(c, t) = \rho c$ dans le cas des glioblastomes [2].

Cette équation modélise deux phénomènes couplés, la diffusion et la multiplication des cellules tumorales dans le parenchyme cérébral. Les plus récentes publications [2] font état de deux constantes de diffusion d différentes pour matière blanche et matière grise, et utilisent une source ponctuelle. Dans cet article, nous proposons deux évolutions de ce modèle : une version anisotrope de l'équation de réaction diffusion et une nouvelle compartimentation du modèle pour tenir compte du centre de la tumeur.

<figure><img src="image_2.png" /><figcaption>FIG. 1 – Évolution de la croissance d'un glioblastome en 6 mois et segmentation utilisée dans les protocoles de radiothérapie. (A) IRM T2 du patient en mars 2001. (B) Segmentation proposée en radiothérapie : le GTV1 (en rouge) correspond à la partie la plus active de la tumeur ; le GTV2 (en bleu) correspond à la zone infiltrante de la tumeur. (C) IRM T2 du même patient en septembre 2001.</figcaption></figure>

### 2.1.1 Compartimentation du modèle en trois zones distinctes

Nous avons isolé pour le glioblastome 3 modes d'envahissement différents dans le parenchyme cérébral :

- La partie centrale du glioblastome, appelé GTV1 (Gross Tumor Volume 1) dans les protocoles de radiothérapie (zone rouge de la figure 1) constitue la partie la plus "active" de la tumeur. La partie périphérique

<footnote n="1">¹"_" représente un vecteur 3 × 1 et "_" une matrice 3 × 3</footnote>
du GTV1 est le lieu de nombreuses multiplications cellulaires. La partie centrale du GTV1 est souvent composée d'un mélange de cellules actives, de cellules mortes et de nécrose. L'ensemble de ce volume est saturé en cellules tumorales, si bien que l'augmentation du nombre de cellules se traduit par une augmentation de son volume.

- La partie périphérique de la tumeur, appelée GTV2 (zone bleue) présente un gradient de cellules tumorales décroissant du centre vers la périphérie. Cette zone est le lieu d'une migration, ou diffusion, des cellules essentiellement à travers les fibres de matière blanche.

- Enfin, la diffusion des cellules du GTV2 est très fortement freinée dans la matière grise, et arrêtée à la frontière du cerveau (crâne et la faux du cerveau -structure inter-hémisphérique- et ventricules).

Notre modèle reprend cette compartimentation : le GTV1 est une zone dont la densité de cellules tumorales est saturée $c = c_{max}$. Les cellules obéissent à la loi de réaction diffusion dans la reste du cerveau, avec des paramètres différents pour matière blanche et matière grise. Cette équation d'évolution de $c$ est également stoppée au delà de $c_{max}$. Enfin les conditions aux limites de Neumann imposent un flux nul à l'interface entre le cerveau et le liquide céphalo-rachidien (LCR) qui l'entoure. La valeur de la constante de densité de saturation est issue de la littérature [1] : $c_{max} = 3.5 \times 10^4 \text{ cellules mm}^{-3}$. Les valeurs des coefficients de multiplications cellulaires $\rho$ et de diffusion dans la matière grise $d$ sont initialisées à partir des données de la littérature, puis ajustées manuellement aux observations chez notre patient : $\rho = 2.2 \times 10^{-5} \text{day}^{-1}$ et $d = 10^{-7} \text{mm}^2 \text{s}^{-1}$.

### 2.1.2 Anisotropie de diffusion

Par définition, le glioblastome se développe en suivant les fibres de matière blanche. Cela implique une vitesse de diffusion différente dans la direction des fibres et dans la direction transverse. Nous proposons d'introduire cette anisotropie à travers un tenseur de diffusion $D_T$ remplaçant le coefficient $d$ dans l'équation 1 et décrivant localement la probabilité de diffusion des cellules tumorales en fonction de la direction (un exemple d'image de diffusion est donné en figure 2, seul le premier vecteur propre du tenseur est représenté). Cette donnée n'étant pas disponible chez le patient, nous faisons l'hypothèse qu'il existe une relation de proportionnalité entre la diffusion des molécules d'eau $D_W$ et la diffusion des cellules tumorales : $D_T = \alpha D_W$. On peut alors utiliser une IRM de diffusion pour initialiser ce champ de tenseurs chez le patient. La valeur de $\alpha$, différente pour chaque tumeur et caractérisant son agressivité, est également ajustée manuellement : $\alpha = 5.10^{-3}$. Cette valeur reste néanmoins cohérente avec les valeurs proposées par Swanson [2].

<figure><img src="image_3.png" /><figcaption>FIG. 2 – IRM de diffusion faisant apparaître les faisceaux de matière blanche.</figcaption></figure>

L'ensemble du modèle d'envahissement des cellules tumorales (compartimentation et anisotropie) est synthétisé sur la figure 3, image de gauche.

## 2.2 Couplage mécanique

Les travaux de Miller [3], uniques dans le domaine de la rhéologie du cerveau, on montré un comportement visco-élastique non linéaire du parenchyme cérébral. On peut cependant négliger les effets non-stationnaires dans notre cas du fait de l'échelle de temps considérée (plusieurs mois). Nous proposons de plus une version linéarisée de cette loi, conduisant à un module de Young $E = 694 \text{ Pa}$ et un coefficient de poisson de 0.45. Nous avons montré [4] que l'erreur ainsi commise sur la contrainte est inférieure à $4 \text{ Pa}$ pour des déformations inférieures à 10 %.

Wasserman a proposé dans [5] de modéliser l'influence mécanique d'une tumeur par une augmentation de pression P proportionnelle à $N/V$, avec N le nombre total de cellules tumorales et V le volume total de la tumeur. Nous proposons ici une nouvelle relation d'équilibre pour le cerveau reliant la contrainte à la densité de cellules tumorales :

$$ \mathrm{div} \left( \underline{\underline{\sigma}} - \lambda \underline{\underline{\Pi}}_3 c \right) + \underline{\underline{f}}_{\mathrm{ext}} = 0 \qquad (2) $$

Cette équation est en réalité la version différentielle de la loi proposée par Wasserman. Elle peut s'interpréter comme une augmentation de la pression locale du milieu proportionnelle à la densité de cellules tumorales infiltrées dans ce milieu. La valeur de la constante $\lambda$ est ajustée de manière automatique en minimisant l'erreur sur le déplacement mesurée à partir des images : $\lambda = 1.4 \times 10^{-9} \text{ N mm Cells}^{-1}$.
Nous appliquons cette équation de comportement dans tout le cerveau, à l'exception de la faux dont les pa-
ramètres de rigidité sont bien plus élevés (200 000 Pa [6]). La surface extérieure du cerveau est fixée contre le
crâne, et les ventricules sont laissés libres de contraintes (cette modélisation repose sur l'hypothèse que la tumeur
n'affecte pas la circulation du LCR de l'espace intra-crânien).

L'ensemble du modèle mécanique est synthétisé sur la figure 3, image de droite.

<figure><img src="image_4.png" /><figcaption>FIG. 3 – Synthèse des modèles de diffusion (gauche) et mécanique (droite)</figcaption></figure>

2.3 Construction du modèle éléments finis à partir d'une image médicale

2.3.1 Construction du maillage

Nous avons pour objectif de pouvoir réaliser une comparaison locale à la fois qualitative et quantitative entre
la simulation issue du modèle proposé et la croissance d'un véritable gliolastome chez un patient donné. Cette
simulation repose sur un maillage tétraédrique représentant la géométrie du cerveau du patient. La construction de
ce maillage est réalisée en plusieurs étapes à partir de l'IRM anatomique du patient :

- La segmentation du cerveau du patient dans l'IRM 3D du patient (256 x 256 x 120 mm), autrement dit la réalisation d'une classification binaire pour chaque voxel (pixel 3D) de l'image du patient en fonction de son appartenance à la classe "cerveau". Ce problème de segmentation du cerveau est désormais un problème classique de traitement d'images, nous avons utilisé pour cela le logiciel *brainvisa* disponible sur internet <fnref n="2" />.

- La triangulation de la surface extérieure (cortex) et intérieure (ventricules) du cerveau. Cette triangulation est extraite à partir de l'iso-surface (continue) de la segmentation précédente à l'aide d'un algorithme issu de l'algorithme de Chew [7] et développé par l'équipe Geometrica de l'INRIA Sophia-Antipolis<fnref n="3" />.

- La tétraédrisation du volume défini par les deux surfaces triangulaires précédentes. Nous utilisons pour cela le logiciel ghs3d, développé par l'équipe gamma de l'INRIA Rocquencourt<fnref n="4" />

Le maillage ainsi généré servira de support au calcul éléments finis. Nous utilisons un maillage relativement fin (300 000 tétraèdres pour un volume d'environ 1400 cm<fnref n="3" />) afin de pouvoir représenter la complexité des sous-structures présentes à l'intérieur du cerveau.

2.3.2 Étiquetage des tétraèdres du maillage

La procédure d'étiquetage des tétraèdres du maillage nécessite une segmentation préalable de l'image du pa-
tient. Nous utilisons ici une méthode automatique reposant sur le recalage d'un atlas sur le patient. Cette méthode
met en correspondance chaque voxel de l'image du patient avec un voxel de l'image de l'atlas portant un label. Ce
label est ensuite reporté sur chaque voxel de l'image du patient. Chaque tétraèdre se voit alors assigné la classe ma-
joritaire des voxels qu'il contient. Nous utilisons une classe de déformations non rigide affine pour recaler l'image
du sujet avec celle de l'atlas. Ce type de transformation à faibles degrés de liberté (12) garantie la robustesse de

<footnote n="2">²www.brainvisa.info</footnote>

<footnote n="3">³www-sop.inria.fr/geometrica</footnote>

<footnote n="4">⁴www-rocq1.inria.fr/gamma/</footnote>
l'algorithme à la présence de tumeur. La contrepartie associée est une correspondance approximative entre la segmentation et l'image du patient. Néanmoins, l'effet de moyenne dû au rapport moyen du nombre de voxels par tétrahèdre (≈ 10) lisse les conséquences d'une segmentation approximative.

## 3 Résultats

Le modèle précédemment décrit est utilisé pour simuler la croissance d'un glioblastome en utilisant les images du patient pour son initialisation. Les résultats sont présentés en deux parties : la diffusion des cellules tumorales et les effets mécaniques.

<figure><img src="image_5.png" /><figcaption>FIG. 4 – 1. IRM T2 03/2001, 2. IRM T2 03/2001 + initialisation du GBM, 3. IRM T2 09/2001, 4. IRM T2 09/2001 + courbes d'iso-densité de cellules de GBM</figcaption></figure>

### 3.1 Diffusion des cellules tumorales

La validation par l'image de la simulation de l'évolution de la densité de cellules tumorales se heurte à deux problèmes :

- L'IRM ne permet pas une observation directe de la tumeur, mais de l'oedème cérébral. Toutefois, la corrélation locale entre oedème et tumeur est largement admise par la communauté médicale.

- L'IRM, par opposition à la tomodensitométrie, n'est pas une technique calibrée, si bien qu'on ne peut mesurer la composition du tissu avec cet outil. Cependant, là aussi, la corrélation spatiale entre position de la tumeur sur l'IRM et le scanner a été vérifiée. Nous utiliserons donc le seuil de détection en tomodensitométrie sur la densité de cellules tumorales (8000 cells mm⁻³) pour la visualisation des courbes d'iso-densité de cellules tumorales.

La figure 4 montre une comparaison entre la simulation et la croissance réelle de la tumeur chez le patient. La première colonne montre l'image du patient en mars 2001 ; la deuxième colonne les contours utilisés pour l'initialisation du modèle. La troisième colonne montre l'image du patient acquise 6 mois plus tard. Elle permet de prendre la mesure de l'ampleur de l'envahissement de la tumeur en 6 mois. La dernière colonne montre les courbes d'iso-densité de cellules tumorales simulées par le modèle, au dessus du seuil de détection $c_{max}$.

Ces résultats qualitatifs montrent une très bonne corrélation entre la simulation et l'observation dans les images. On retrouve en outre les "branches" d'infiltration de la tumeur autour des noyaux gris centraux.

### 3.2 Effet mécanique

L'image 5 présente les déplacements simulés par le modèle. Le déplacement simulé moyen de l'interface GTV1-GTV2 est de 3mm. On observe également un effet classique de compression des ventricules avec une
variation du volume des ventricules latéraux de 4.6 ml. Afin de quantifier l'erreur sur le déplacement, un expert médical a manuellement sélectionné des points homologues dans les images de mars et septembre. Le déplacement ainsi mesuré est comparé à celui simulé par le modèle. L'erreur moyenne sur le déplacement est de 1.3 mm pour un déplacement moyen sur les points de 2.7 mm. Sans être capable de capturer l'intégralité des déplacements réels, cette simulation démontre qu'une telle modélisation est capable d'expliquer une grande partie du déplacement. L'erreur résiduelle peut être due à différents facteurs :

- Le rapport entre le déplacement moyen (2.7 mm) et la taille de chaque voxel (1 mm) est dans la marge d'erreur liée à la sélection manuelle de point homologues.

- Le maillage réalisé ne comprend pas l'espace intersticiel des sillons qui peut subir une compression plus grande.

<figure><img src="image_6.png" /><figcaption>FIG. 5 – Déplacement du tissus induit par l'effet de masse de la tumeur. (Gauche) Vue subjective 3D des déplacement sur le modèle. (Droite) projection des déplacements sur une coupe de l'image.</figcaption></figure>

## 4 Conclusion

Nous avons proposé dans cet article plusieurs évolutions des modèles de croissance de tumeurs cérébrales existants, à la fois dans l'équation d'évolution de la densité de cellules tumorales, dans la compartimentation spatiale des phénomènes, et dans le couplage physio-mécanique. À ces contributions s'ajoute la comparaison quantitative avec les images d'un patient.

Les premiers résultats présentés ici sont très encourageants, aussi bien du point de vue envahissement que mécanique. Ils ouvrent désormais une voie vers une meilleure quantification des paramètres, et notamment vers l'étude des paramètres spécifiques au patient et/ou à la tumeur. Des études portant sur un grand nombre de sujet sont néanmoins nécessaires à la validation d'un tel modèle. Les nouvelles techniques d'imagerie telle que l'imagerie spectroscopique constituent certainement un axe de recherche privilégié pour une telle validation in-vivo non-invasive.

## Références

[1] P. Tracqui, G. Cruywagen, D. Woodward, G. Bartoo, J. Murray, and E. Alvord Jr, "A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth," *Cell Proliferation*, vol. 28, no. 1, pp. 17-31, Jan 1995.

[2] K. Swanson, E. Alvord Jr, and J. Murray, "Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy," *British Journal of Cancer*, vol. 86, no. 1, pp. 14-18, Jan 2002.

[3] K. Miller, *Biomechanics of Brain for Computer Integrated Surgery*. Warsaw University of Technology Publishing House, 2002.

[4] O. Clatz, P. Bondiau, H. Delingette, M. Sermesant, S. Warfield, G. Malandain, and N. Ayache, "Brain tumor growth simulation," INRIA, Research report 5187, 2004. [Online]. Available : http://www-sop.inria.fr/rapports/sophia/RR-5187.html

[5] R. Wasserman and R. Acharya, "A patient-specific in vivo tumor model," *Mathematical Biosciences*, vol. 136, no. 2, pp. 111-140, Sep 1996.

[6] M. Schill, M. Schinkmann, H.-J. Bender, and R. Männer, "Biomechanical simulation of the falx cerebri using the finite element method," in *18. Annual International Conference, IEEE Engineering in Medicine and Biology. S*, 1996.

[7] L. P. Chew, "Guaranteed-quality mesh generation for curved surfaces," in *SCG '93 : Proceedings of the ninth annual symposium on Computational geometry*. ACM Press, 1993, pp. 274-280.