\documentclass[XUPS,XML,SOM,Unicode,francais, NoFloatCountersInSection, NoEqCountersInSection]{cedram}
\usepackage{xups22-02}
\setcounter{tocdepth}{2}
\let\subsection\Subsection
%\XUPScorrections

\hyphenation{notion noter résis-tant résis-tants modè-le modè-les tumeur tumeurs}

\begin{document}
\frontmatter

\title[Modélisation mathématique des gliomes de bas grade]{Modélisation mathématique des\nobreakspace gliomes\nobreakspace de\nobreakspace bas\nobreakspace grade}

\author[\initial{E.} \lastname{Grenier}]{\firstname{Emmanuel} \lastname{Grenier}}
\address{Unité de Mathématiques Pures et Appliquées et CNRS, UMR 5669, École Normale Supérieure de Lyon,
46, allée d’Italie, 69364 Lyon Cedex 07, France}
\email{emmanuel.grenier@ens-lyon.fr}
\urladdr{http://perso.ens-lyon.fr/egrenier/}

\author[\initial{B.} \lastname{Ribba}]{\firstname{Benjamin} \lastname{Ribba}}
\address{Roche, Bâle, Suisse}
\email{}
\urladdr{}

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


\begin{abstract}
Ce texte est une introduction à la modélisation mathé\-matique en cancérologie, et en particulier à la modélisation des gliomes, qui sont une forme particulière de tumeurs cérébrales. L'objec\-tif est d'illustrer comment une modélisation mathématique peut contribuer à répondre à des questions médicales et à mieux comprendre comment traiter des tumeurs cancéreuses.
\end{abstract}
\maketitle
\tableofcontents

\mainmatter
\section{Introduction}

\subsection{Gliomes}

L'objectif de cette section est de présenter rapidement et de façon très simplifiée ce qu'est un gliome.
La population des patients atteints de gliomes est très hétérogène, et l'évolution de ces tumeurs varie très fortement
d'un patient à l'autre. En particulier les chiffres donnés ci\nobreakdash-dessous ne sont que des moyennes, données
pour une sous catégorie de gliomes.

Les cellules gliales sont des cellules du cerveau qui servent de \og support\fg aux neurones en leur apportant des nutriments
et de l'\hbox{oxygène}, en éliminant les déchets et en produisant la myéline qui entoure les fibres nerveuses.
Suite à diverses mutations, ces cellules gliales peuvent se multiplier de façon incontrôlée, conduisant à des
gliomes, ou tumeurs gliales.

Les gliomes représentent $80 \%$ des tumeurs cérébrales malignes. Ils sont généralement détectés après des
crises d'épilepsie ou des signes de pression intra-crânienne trop élevée. Le diagnostic fait suite à des examens de type
IRM, voire à des biopsies. Les gliomes sont classifiés en \og grades\fg, du grade $1$ au grade $4$, du moins sévère
au plus sévère. Le grade $1$ correspond
à des gliomes curables par chirurgie, le grade $2$ à des gliomes de bas grade, diffus et infiltrants. Les grades $3$
et $4$ (haut grades) sont plus agressifs, avec hélas des durées de survie moindres.
Les gliomes sont de plus caractérisés par leur profil génétique, c'est-à-dire en fonction des mutations ayant eu lieu
(co-délétion $1$p/$19$q, mutation TP$53$, mutation IDH,...).

Cet exposé porte sur la modélisation des gliomes de bas grades. Pour ces types de gliomes, les cellules restent bien
différenciées: les cellules tumorales ont les caractéristiques des cellules cérébrales. Elles se divisent peu, toutefois
elles sont en excès et forment une masse anormale dans le cerveau. L'évolution est lente et peut s'étendre sur~$5$ voire $10$ ans, mais leur évolution vers des tumeurs de grade $3$ ou $4$ est inéluctable.

Les traitements se divisent en trois catégories: la chirurgie, la radiothérapie et la chimiothérapie.
La chirurgie qui consiste à tenter d'enlever la tumeur se heurte à deux difficultés: d'une part la tumeur n'a pas de bords
précis, car les cellules tumorales \og diffusent\fg dans le tissu environnant. Il en reste toujours après chirurgie et la rechute
est inéluctable. D'autre part, certaines zones du cerveau ne peuvent pas être ôtées sans compromettre la qualité de vie
du patient.
La radiothérapie est un traitement standard, parfois en combinaison avec une chirurgie.

\begin{figure}[htb]
\centerline{\includegraphics[width=11cm]{xups22-02_figures/figure1.png}}
\caption{\small
Évolution d'un gliome de bas grade traité par PCV. Le graphique de gauche représente le diamètre moyen au cours du temps
(années), le traitement de PCV est représenté par la bande noire \og CT\fg. À droite, IRM de ce patient aux temps E, F, G et H. }
\label{figure1}
\end{figure}

\smallskip

La chimiothérapie est de plus en plus utilisée car elle induit des répon\-ses tumorales efficaces et peut stabiliser
la maladie. Elle peut être donnée en premier traitement, peut permettre de libérer des zones cérébrales sensibles
avant chirurgie ou peut permettre de repous\-ser l'utilisation de la radiothérapie.

\smallskip

Dans cet exposé nous nous intéressons à la modélisation de la chimiothérapie.
Il existe deux chimiothérapies différentes: le PCV (procarbaine, lomustine et vincristine)
et le témozolomide (TMZ).

Le PCV est une combinaison de trois molécules, administrées en six cycles maximum car il est très toxique.
Chaque cycle dure six semai\-nes, sous divisé en administration des diverses molécules aux jours $1$, $8$, $21$ et $29$.
La figure \ref{figure1} montre l'évolution du diamètre moyen de la tumeur chez un patient. La maladie est détectée en E et le
patient reste sans traitement pendant $6$ ans. Pendant ce temps le diamètre passe de $35$ à $60$ mm. Puis
lors de la chimiothérapie ce diamètre tombe à $45$ mm avant de poursuivre sa décroissance pendant presque
trois ans (jusqu'à atteindre $35$ mm)
puis de recroître à nouveau (\hbox{rechute}). Cette décroissance prolongée est particulièrement étonnante
car les molécules actives sont très rapidement éliminées par l'organisme, et~leur concentration est nulle après quelques
jours.

Le TMZ est moins toxique que le PCV et les patients peuvent être traités plus longtemps avec cette molécule (jusqu'à
douze cycles).\enlargethispage{1.5\baselineskip}%

\subsection{Données cliniques}

Les données cliniques sont des données d'imagerie (IRM). Chaque patient passe régulièrement un IRM,
soit pour suivre la croissance avant traitement, soit pour suivre l'efficacité du traitement, soit pour suivre la rechute.
Les données présentées proviennent des hôpitaux de Lyon (Hôpital Pierre Wertheimer
et Centre Léon Bérard) et concernent $21$ patients pour le PCV, $109$ pour le TMZ et $12$ pour la radiothérapie.

L'IRM fournit des images 3D du cerveau du patient que le clinicien peut analyser pour localiser et mesurer la tumeur.
La tumeur est également tridimensionnelle mais il est compliqué et long d'en tracer les contours.
Comme la tumeur a une forme généralement ovoïde, les cliniciens se contentent dans leur pratique quotidienne
de mesurer deux \og diamètres\fg (le plus grand et le plus petit) et de les combiner
pour en tirer un taille tumorale (en mm). Le volume est aussi évalué à partir de ces deux diamètres.

Ce n'est bien sûr guère précis, mais c'est rapide, et d'autre part l'interprétation précise des IRM est délicate.
Deux cliniciens différents ne vont pas exactement \og voir\fg la tumeur au même endroit car ses contours sont flous.
Le même clinicien à qui
le même IRM est présenté deux fois de suite ne va pas contourer la tumeur exactement au même endroit.
Certaines zones de l'IRM sont ambiguës. L'écart de volume mesuré entre deux cliniciens peut atteindre $20\%$.
Par ailleurs, pour donner une interprétation clinique de l'évolution de la taille d'une tumeur
les cliniciens utilisent le critère dit \og RECIST\fg, qui indique
que l'on peut dire qu'une tumeur croît si son volume augmente d'au moins $20\%$ et qu'elle décroît si son volume
décroît d'au moins $30\%$. Entre $-30\%$ et $+20\%$ la tumeur est dite sans évolution ou stable.
Ce~critère est cohérent avec la faible précision de mesure du diamètre tumoral.

Il est très important de souligner que les données cliniques sont toujours entachées d'erreurs de mesure considérables,
ce qui contraste fortement avec la physique ou la chimie où il n'est pas rare d'avoir des marges d'erreurs inférieures au pourcent,
voire beaucoup moins, ce qui conduit à la notion de \og variabilité intrapatient\fg.\enlargethispage{-2.5\baselineskip}%

Les figures \ref{figure3}, \ref{figure4} et \ref{figure5} montrent diverses données cliniques.
La principale caractéristique de ces données est l'extrême variabilité d'un patient à l'autre, appelée
\og variabilité interpatient\fg. Chez certains patients
l'évolution est lente, chez d'autres elle est nettement plus rapide. Par exemple, sur la figure \ref{figure5},
le patient $2$ répond bien au TMZ, alors que cette même molécule ne fait que freiner l'évolution du gliome chez le patient
$14$. Le patient $47$ rechute même pendant la chimiothérapie. L'évolution d'une tumeur dépend en effet de ses caractéristiques
génétiques (mutations présentes), de l'existence de clones résistants au TMZ avant le début du traitement, ou
de leur apparition lors du
traitement, mais aussi de la tolérance de la chimiothérapie par le patient (effets secondaires) et de l'état général du patient.
Contrairement à la physique ou à la chimie, des données obtenues dans des circonstances similaires peuvent être fortement
différentes car elles dépendent du patient étudié.
Ceci conduit à la notion fondamentale de variabilité interpatient (d'un patient à l'autre), notion que nous retrouverons
plus loin.

\begin{figure}[htb]
\centerline{\includegraphics[width=10cm]{xups22-02_figures/figure3.png}}
\caption{\small
Dynamique tumorale des $21$ patients traités par PCV. La partie bleue correspond à la croissance naturelle de la
tumeur, la partie verte à la décroissance sous traitement, la partie rouge à la décroissance après traitement et la partie noire
à la recroissance}
\label{figure3}
\end{figure}

\begin{figure}[htb]
\centerline{\includegraphics[width=10cm]{xups22-02_figures/figure4.png}}
\caption{\small
Dynamique tumorale des $77$ patients traités par TMZ. Les courbes rouges correspondent au percentiles
$5$, $50$ et $95$.}
\label{figure4}
\end{figure}

\begin{figure}[htb]
\centerline{\includegraphics[width=9cm]{xups22-02_figures/figure5.png}}
\caption{\small
Illustration des différentes dynamiques observées dans la population traitée par TMZ}
\label{figure5}
\vspace*{-\baselineskip}
\end{figure}

\section{Quelques modèles de dynamique des populations}

Cette partie est consacrée à la présentation de quelques modèles simples et classiques de pharmacométrie
et de croissance tumorale.
La pharmacométrie est la science de l'étude de l'interaction entre l'organisme et un médicament. Elle se décompose en
pharmacocinétique (PK: étude de l'évolution du médicament dans l'organisme) et pharmacodynamie (PD: étude
de l'impact du médicament sur l'organisme, en termes d'efficacité ou de toxicité). La pharmacocinétique - pharmacodynamie
(PK - PD) décrit l'évolution de la réponse au traitement au cours du temps.\enlargethispage{2\baselineskip}

\subsection{Modélisation PK}

Les modèles de pharmacocinétique permettent d'évaluer l'évolution d'un médicament dans l'organisme.
Ce domaine est très étudié par les grandes compagnies pharmaceutiques où des services entiers y sont consacrés.
La PK doit décrire l'absorption, la distribution, la métabolisation et l'élimination des molécules actives. Elle fournit
des données cruciales telles que la demi-vie d'une molécule dans l'organisme.

L'organisme est généralement divisé en \og compartiments\fg où la concentration du médicament est supposée uniforme,
et ne dépendre que du temps (par exemple: le sang, un organe cible, la concentration dans le foie). Les modèles de PK
sont des systèmes d'équations différentielles ordinaires, plus ou moins complexes, décrivant le passage des molécules
d'un compartiment à l'autre et leur dégradation.

Le modèle PK standard du TMZ est complexe, donc nous ne présenterons ici qu'une variante très simplifiée.
Ce modèle simplifié a deux compartiments, le sang et l'organe cible (ici le cerveau), et deux concentrations $G$ (dans le sang)
et $A_1$ (dans le cerveau). L'évolution de $G$ et de $A_1$ est simplement donné par les équations linéaires suivantes:
\begin{align}
\label{PK1}
\frac{dG}{dt} &= - k_a G + \phi(t),
\\
\label{PK2}
\frac{dA_1}{dt} &= k_a G - k_e A_1,
\end{align}
où $\phi(t)$ est le débit d'injection (quantité de TMZ injectée par unité de temps).
Ces équations peuvent se lire de la façon suivante: le TMZ injecté entre dans le compartiment sanguin.
Un flux $k_a G$ de TMZ passe du sang au cerveau par unité de temps. Ce flux est proportionnel à la concentration
de TMZ dans le sang, avec une constante de proportionnalité $k_a$. Ce flux se retrouve dans \eqref{PK2}.
D'autre part le TMZ se dégrade dans le cerveau (terme $- k_e A_1$) avec une demi vie $\ln(2)/k_e$ (voir figure \ref{figure2}).\enlargethispage{\baselineskip}%

\begin{figure}[htb]
\centerline{\includegraphics[width=9cm]{xups22-02_figures/figure2.png}}
\caption{\small
Modèle PK pour le TMZ, évolution de la concentration $A1$ au cours du temps suite à un bolus (dose administrée
rapidement) }
\label{figure2}
\end{figure}

\subsection{Modélisation PD}

Il s'agit ici de modéliser l'effet d'une molécule en en connaissant la concentration. Le modèle le plus simple est appelé
\og loi du $E_{\max}$\fg ou loi de Hill et exprime l'effet $E$ correspondant à une concentration~$C$ de la molécule par la relation\vspace*{-3pt}
\beq \label{Emax}
E = E_{\max}\, \frac{C}{EC_{50} + C}.
\eeq
Dans cette équation, $E_{\max}$ est l'effet maximum et $EC_{50}$ la concentration conduisant à un effet moitié
du maximum. Cette loi est d'une utilisation extrêmement fréquente en PD.

\subsection{Modèles simples de croissance tumorale}

Pour modéliser la croissance des gliomes il faut non seulement un modèle PK-PD du traitement, mais aussi un modèle
de croissance tumorale. Commençons par quelques modèles simples de croissance tumorale.
Soit $y(t)$ la taille de la tumeur (diamètre ou volume). Le~modèle le plus simple est bien sûr\vspace*{-3pt}
\beq \label{Euler}
\dot y = \lambda y
\eeq
qui traduit le fait que les cellules de la tumeur se divisent régulièrement, et qui conduit à une croissance exponentielle
de la tumeur
\[
y(t) = y(0) e^{\lambda t}.
\]
Le temps de doublement de la tumeur est $\ln(2) / \lambda$. Ce modèle n'est toutefois pas réaliste car heureusement
la croissance de tumeurs n'est pas aussi rapide. Les cellules tumorales lorsqu'elles sont nombreuses se heurtent à des
difficultés d'approvisionnement en oxygène et en nutri\-ments. La mortalité des cellules tumorales augmente avec
la taille de la tumeur, ce qui conduit à raffiner le modèle \eqref{Euler} en
\begin{equation}
\label{Gomp}
\begin{aligned}
\dot y &=\lambda y,
\\
\dot \lambda &= - \alpha \lambda.
\end{aligned}
\end{equation}
Dans ce modèle, la vitesse de croissance $\lambda$ décroît au cours du temps pour traduire la difficulté d'accès aux
ressources. La solution de~\eqref{Gomp} est explicite et porte le nom de loi de Gompertz:
\beq \label{Gomp3}
y(t) = y(0) \exp \bigl( C_0 (1 - e^{-C_1 t}) \bigr)
\eeq
où $C_0$ et $C_1$ sont deux constantes. Cette loi est un grand classique de la cancérologie.

Une autre approche est de remarquer que la tumeur ne peut pas dépasser une taille maximale.
Supposons pour simplifier qu'elle soit sphérique, de rayon $R(t)$. Alors ses besoins énergétiques
évoluent comme $R^3(t)$ alors que les apports d'énergie (oxygène et nutriments) qui sont apportés par l'extérieur
doivent transiter par sa surface, et~donc évoluent comme $R^2(t)$. Le rayon ne peut donc pas croître indéfiniment,
et ne peut dépasser une valeur seuil $R_0$ (d'autre part bien sûr un gliome ne peut pas avoir une taille supérieure à la
boite crânienne). Une façon simple de traduire cette limitation est d'écrire
\beq \label{logistique}
\dot R= \lambda R( 1 - \sfrac{R}{R_0}) ,
\eeq
où $\lambda$ est la vitesse de croissance maximale et $R_0$ le rayon maximal.
Ce~modèle porte le nom de modèle logistique et a été introduit dans les années 1830 par Verhulst
pour tenir compte de la limitation naturelle de taille dans des modèles d'écologie.
Il est d'un usage constant en écologie (voir en particulier \cite{Murray}).

À noter que le modèle de Gompertz peut se mettre sous la forme voisine
\beq \label{Gomp4}
\dot y = - \alpha y \log(\theta^{-1} y),
\eeq
où $\theta$ est la taille maximale que la tumeur peut atteindre.

\subsection{Remarques}

Il est important de noter qu'il n'existe pas un seul modèle de croissance tumorale. Contrairement à la physique
ou à la chimie, il n'existe pas de grands principes (lois de la thermodynamique, loi de Newton, principes de la
mécanique quantique, équations de la relativité ou encore équations de Maxwell) à partir desquelles il serait possible
d'étudier un système biologique (quitte à ensuite faire des approximations ou des simplifications).
Les systèmes biologiques sont complexes à cause du grand nombre d'agents en jeux: multiples types de cellules
différentes (cellules saines, tumorales, système immunitaire), multiples échelles spatiales (organe, cellule, noyau cellulaire)
et temporelles (temps de croissance d'une cellule / temps d'apparition de mutations).
Le modèle utilisé dépend de la question étudiée, en suivant un principe de parcimonie: il est toujours préférable
d'étudier le modèle le plus simple possible. En particulier pour se confronter à des données cliniques, le modèle
ne doit pas comporter trop de paramètres sinon ceux-ci ne peuvent pas être évalués correctement car il y a trop de possibilités
de compensations entre leurs divers effets.

\smallskip

D'autre part une tumeur n'est pas homogène. Auparavant les biologistes pensaient qu'une tumeur était relativement homogène,
toutes ses cellules étant similaires, issues d'une même cellule tumorale primitive. On sait, depuis, que les tumeurs sont hautement
hétérogènes et composées de multiples clones, en perpétuelle évolution: apparition de nouveaux clones suite à des mutations,
disparition des clones les moins adaptés, interaction entre les divers clones (il peut arriver que des cellules tumorales moins
agressives \og calment\fg des cellules plus agressives). Le paradigme actuel est de regarder une tumeur comme un écosystème
en évolution (au sens de l'évolution darwinienne),
qu'il faut contrôler le plus longtemps possible en intervenant intelligemment dessus.
Par exemple, un traitement agressif peut être contre-indiqué car il peut éliminer les clones peu agressifs et favoriser ainsi
les clones agressifs qui ont plus de place et accès à plus de ressources,
conduisant à une rechute plus brutale et plus rapide qu'avec un traitement moins agressif.

\section{Un premier modèle d'évolution de gliome}

Ce travail fait suite à une question de F.\,Ducray, cancérologue
lyonnais qui cherchait à comprendre pourquoi la tumeur continuait à décroî\-tre \emph{après} l'arrêt du traitement, ce qui est
inhabituel. La tumeur décroît alors qu'il n'y a plus de tout de substance active dans l'organisme.
F.\,Ducray a émis l'hypothèse de l'existence de deux types de cellules tumorales: des cellules proliférantes
qui, comme leur nom l'indique, se multiplient rapidement, et des cellules quiescentes (qui \og dorment\fg, et, bien que tumorales,
ne se multiplient pas).
La question posée est alors simple:

\smallskip

\noindent \emph{Est-ce que ce scénario est compatible avec les données cliniques?
\hbox{Est-ce} qu'il permet d'expliquer cette décroissance prolongée?}

\smallskip

\noindent Pour cela, un modèle mathématique
de gliome a été développé pour quantifier ce scénario et le comparer aux données patients.

\subsection{Modélisation}

Nous présentons ici le modèle développé dans \cite{Ribba}. Ce modèle considère trois types de cellules: des cellules
proliférantes, des cellules quiescentes non endommagées et des cellules quiescentes endommagées.
Nous noterons $P$, $Q$ et $Q_P$ leurs nombres respectifs. Le schéma~\ref{figure6} donne une
vue synthétique des liens entre ces types cellulaires.

\begin{itemize}

\item Les cellules proliférantes se multiplient avec une vitesse $\lambda_P$, meurent sous l'effet du traitement avec une vitesse
maximale $\gamma$, et~une certaine proportion devient quiescente ($k_{PQ}P$).

\item Les cellules quiescentes ne se multiplient pas mais sont endommagées sous l'effet du traitement.

\item Les cellules quiescentes endommagées meurent ou sont réparées (mutations). Dans ce second cas, elles deviennent
à nouveau proliférantes (ce qui modélise l'apparition de résistance).

\item Il faut aussi modéliser la concentration du PCV et son effet. Pour simplifier on suppose que le PCV a une demi-vie
$\ln(2)/k_{DE}$ et un effet proportionnel à sa concentration. On note $\phi(t)$ la quantité de PCV injectée par unité de
temps.

\end{itemize}

Ceci conduit au système suivant
\begin{align*}
\frac{dC}{dt} &= -k_{DE} C + \phi(t)
\\
\frac{dP}{dt}&= \lambda_P P \Bigl( 1- \sfrac{P^\star}{K} \Bigr) + k_{Q_PP} Q_P - k_{PQ} P -\gamma_P k_{DE} C P
\\
\frac{dQ}{dt} &= k_{PQ}P- \gamma_Q k_{DE} C Q
\\
\frac{dQ_P}{dt} &=\gamma_Q k_{DE} C Q - k_{Q_P P} Q_P -\delta_{QP} Q_P
\\
P^\star &= P + Q + Q_P.
\end{align*}
À noter que $P^\star$ est la taille totale de la tumeur, et que $K$, appelée capacité totale du milieu, est la taille
maximale possible de la tumeur.

\begin{figure}[htb]
\centerline{\includegraphics[width=10cm]{xups22-02_figures/figure6.png}}
\caption{Schéma du modèle de gliome }
\label{figure6}
\end{figure}

\subsection{Paramétrisation}

Ce système comporte de très nombreux paramètres: $k_{DE}$ (\hbox{vitesse} de dégradation du PCV),
$\lambda_P$ (vitesse maximale de croissance du gliome), $K$ (taille maximale du gliome), $k_{Q_P P}$ (réparation),
$k_{PQ}$ (transition de $P$ à $Q$), $\gamma_P$ (efficacité du PCV sur $P$), $\gamma_Q$ (efficacité sur $Q$),
$\delta_{QP}$ (mortalité de $Q_P$), soit $8$ paramètres au total.
À ceux-ci doivent s'ajouter les conditions initiales, à savoir $P$, $Q$, $Q_P$ lors du premier IRM, ce qui fait $11$
paramètres. Nous rassemblons ces paramètres en un vecteur $\theta$ à $11$ composantes.

Si $\theta$ est donné, alors il est possible de calculer une solution appro\-chée du système ci-dessus pour tout temps $t$,
en utilisant une méthode numérique d'intégration des équations différentielles ordinaires (\hbox{Euler} ou Runge-Kutta par exemple).
Soit $F(t,\theta)$ la valeur de la taille totale~$P^\star$ de la tumeur au temps $t$ avec les paramètres $\theta$.

\Changel

Une première approche est de chercher à ajuster ces paramètres patient par patient. Plus précisément, pour
chaque patient $P_i$ nous disposons de mesures de diamètres $d_{i,l}$ pris à des temps $T_{i,l}$.
Une idée naturelle est de chercher le jeu de paramètres $\theta$ qui rend compte au mieux des observations, c'est-à-dire
de trouver $\theta$ qui minimise l'écart quadratique $J_i$ entre les observations et les prédictions du modèle:
\beq \label{defiJ}
J_i(\theta) = \sum_l \bigl| F(T_{i,l},\theta) - d_{i,l} \bigr|^2.
\eeq
Minimiser $J_i$ est une approche d'ajustement de paramètres (aussi appe\-lée régression non linéaire ou \og nonlinear fitting\fg
en anglais). Le paramètre $\theta_i$ qui minimise $J_i$ est obtenu numériquement en utilisant un algorithme
de minimisation de fonction, par exemple un algorithme de gradient à pas fixe.

Toutefois cette approche échoue car il y a généralement moins de données par patient que de paramètres.
Il faudrait au moins $11$ IRM par patient pour que cette approche ait un sens, et encore, au vu des marges d'erreur,
même avec $11$ IRM l'ajustement de paramètres \og \hbox{patient} par patient\fg a peu de sens.
En général il y a moins de mesures que de paramètres et ces derniers ne peuvent pas être identifiés.

\subsection{Approche de population}

La seule possibilité de paramétrer ce modèle est de considérer toutes les données de tous les patients
simultanément, au lieu de considérer chaque patient individuellement. Cette approche est appe\-lée \og approche
populationnelle\fg.

\smallskip

La première étape est de décrire l'erreur de mesure, c'est-à-dire de développer un \og modèle d'erreur\fg.
Imaginons que le vrai diamètre de la tumeur soit $d_0$. Un IRM interprété par un clinicien va donner une mesure
$d_\mathrm{mes}$, a priori différente de $d_0$. La différence $d_\mathrm{mes} - d_0$ est l'erreur de mesure.
La première étape consiste à faire une hypothèse sur cet écart, c'est-à-dire à décrire la loi de probabilité qu'il suit.
Par exemple, il est souvent supposé que l'erreur est gaussienne, de~moyenne nulle et d'écart type $\sigma$
\beq \label{erreur}
d_\mathrm{mes} - d_0 \sim {\cal N}(0,\sigma),
\eeq
autrement dit la probabilité $P(d_\mathrm{mes} - d_0)$ de faire une erreur de mesure $d_\mathrm{mes} - d_0$ est
\beq \label{erreur1}
P(d_\mathrm{mes} - d_0) = \frac{1}{\sigma \sqrt{2 \pi}}\, e^{- (d_\mathrm{mes} - d_0)^2 / 2 \sigma^2}.
\eeq
Il existe bien sûr des modèles d'erreur plus compliqués. L'évaluation de $\sigma$ est complexe car $d_\mathrm{mes}$
est inconnue.

\smallskip

La seconde étape est de prendre en compte la variabilité interpatient (c'est-à-dire la variabilité
de comportement entre un patient et un autre).
Pour cela nous supposons que chaque paramètre suit une loi de probabilité dans la population.
Par exemple on suppose que le paramètre $k_{DE}$ est réparti dans la population suivant une loi gaussienne
de moyenne $k_{DE}^m$ et d'écart type $k_{DE}^\sigma$
\beq \label{kDE}
k_{DE} \sim {\cal N}(k_{DE}^m,k_{DE}^\sigma).
\eeq
De telles hypothèses sont faites sur toutes les composantes de $\theta$
\beq \label{theta}
\theta \sim {\cal N}(\theta^m,\theta^\sigma)
\eeq
où $\theta^m$ et $\theta^\sigma$ sont des vecteurs de $11$ composantes. À noter qu'il peut y avoir des covariances
entre ces différents paramètres, ce que nous ne supposerons pas ici.

Au lieu de chercher $11$ paramètres par patient, nous sommes amenés à chercher $2 \times 11 +1=23$ paramètres globaux
($2$ par composante de $\theta$ plus la valeur de $\sigma$).
Au lieu d'avoir une petite dizaine de données par patient, conduisant à une
paramétrisation mal posée, nous disposons de toutes les données patients (plusieurs centaines) pour $23$ paramètres,
ce qui conduit à un problème de paramétrisation bien conditionné.

Reste à décrire comment chercher ces paramètres \og de population\fg.
Imaginons que $\theta^m$ et $\theta^\sigma$ sont connus. Considérons le patient~$i$. Alors la probabilité que
ce patient ait un jeu de paramètres $\theta_i= (\theta_{i,j})_{1 \le j \le 11}$ est
\[
\prod_{j=1}^{11} \frac{1}{\theta_j^\sigma \sqrt{2 \pi} }\, e^{- ( \theta_{i,j} - \theta_j^m)^2 / 2 (\theta_j^\sigma)^2}.
\]
Si le paramètre du patient est $\theta_i$ alors la probabilité d'observer les données $d_{i,j}$ est
\[
P \bigl( (d_{i,l})_l \mid \theta_i \bigr)\! =\! \prod_l\! \frac{1}{\sigma \sqrt{2 \pi}}\, e^{ - (d_{i,l} - F(T_{i,l},\theta_i))^2 / \sigma^2}
\!=\! \Bigl( \frac{1}{\sigma \sqrt{2 \pi}} \Bigr)^{N_i}\! e^{- J_i(\theta_i) /\sigma^2},
\]
où $N_i$ est le nombre de données cliniques pour le patient $i$.
À~noter que $J_i$, défini en \eqref{defiJ} apparaît dans l'exponentielle. L'approche du para\-graphe précédent
est de minimiser $J_i$, ce qui revient à maximiser la probabilité $P( (d_{i,l})_l \mid \theta_i)$
d'observer les données sachant les paramètres, sans prendre en compte la population des autres patients.

La probabilité d'observer le jeu de données $(d_{i,l})_l$ chez le patient $i$ en prenant en compte la statistique
de $\theta$ dans la population est le produit des deux probabilités précédentes, à savoir
\begin{multline}
\label{patient}
\int \prod_{j=1}^{11} \frac{1}{\theta_j^\sigma \sqrt{2 \pi} }\, e^{- ( \theta_{i,j} - \theta_j^m)^2 / 2 (\theta_j^\sigma)^2}\\[-12pt]
{}\cdot
\prod_l \frac{1}{\sigma \sqrt{2 \pi}}\, e^{ - (d_{i,l} - F(T_{i,l},\theta_i))^2 / \sigma^2} \, \, d\theta_i.
\end{multline}
Les patients étant supposés indépendants, la probabilité d'observer le jeu de données $D = (d_{i,l})_{i,l}$
connaissant la valeur de $\theta^m$ et $\theta^\sigma$ est
\begin{multline}
\label{PD}
P(D \mid \theta^m,\theta^\sigma) =\prod_i
\int \prod_{j=1}^{11} \frac{1}{\theta_j^\sigma \sqrt{2 \pi} }\, e^{- ( \theta_{i,j} - \theta_j^m)^2 / 2 (\theta_j^\sigma)^2}\\[-5pt]
{}\cdot
\prod_l \frac{1}{\sigma \sqrt{2 \pi}}\, e^{ - (d_{i,l} - F(T_{i,l},\theta_i))^2 / \sigma^2} \, \, d\theta_i.
\end{multline}
Cette expression donne la probabilité d'observer l'ensemble des données cliniques $D$ de tous les patients
en fonction de $\theta^m$ et de $\theta^\sigma$, or ce qui nous intéresse
est la probabilité d'avoir $\theta^m$ et $\theta^\sigma$ comme paramètres de population, connaissant les mesures $D$.
Il faut donc \og renverser\fg le point de vue et évaluer $P(\theta^m,\theta^\sigma \mid D)$. Pour cela on utilise la loi de Bayes
\[
P(A |B) P(B) =P(B| A) P(A) =P(A \cap B).
\]
La probabilité d'avoir $B$ sachant $A$ multipliée par la probabilité d'avoir $A$ est égale à la probabilité d'avoir
$A$ sachant $B$ multipliée par $B$, et vaut la probabilité d'avoir simultanément $A$ et $B$.

En utilisant cette formule il vient
\beq\label{Bayes}
P( D \mid \theta^m,\theta^\sigma) P(\theta^m,\theta^\sigma) =
P(\theta^m,\theta^\sigma \mid D) P(D)
\eeq
Le terme $P(\theta^m,\theta^\sigma)$ décrit la probabilité d'avoir $\theta^m$ et $\theta^\sigma$ \og dans l'absolu\fg, c'est-à-dire
avant de connaître les données. Il porte le nom de \og priori\fg dans le vocabulaire de la statistique bayesienne.
Ici, comme rien n'est connu a priori, $P(\theta^m,\theta^\sigma) = 1$ (qui n'est pas vraiment une probabilité...).
Le terme $P(D)$ est la probabilité d'observer ce jeu de données. Il est inconnu, mais ce n'est pas gênant car il peut
être vu comme la constante qu'il faut pour que $P( D \mid \theta^m,\theta^\sigma)$ soit d'intégrale~$1$.
Nous avons donc
\beq \label{Bayes2}
P(\theta^m,\theta^\sigma \mid D) = Z^{-1} P( D \mid \theta^m,\theta^\sigma)
\eeq
où $Z$ (qui vaut $P(D)$) est juste une constante de normalisation.

Il faut maintenant chercher les paramètres de population $\theta^m$ et $\theta^\sigma$ qui minimisent cette expression complexe.
Heureusement des algorithmes efficaces existent, en particulier l'algorithme dit \og SAEM\fg \cite{SAEM} (simulated annealing expectation
maximisation) qui combine de superbes idées venant des chaînes de Markov et de l'algorithme de Metropolis Hastings.
Cet algorithme fournit les paramètres de popu\-lation les plus probables en quelques minutes sur un ordinateur portable.
À noter que \eqref{Bayes2} fournit en fait une mesure de probabilité sur~$\theta^m$ et~$\theta^\sigma$, mesure
qui peut être réduite à son mode le plus probable pour simplifier, mais qui peut être considérée dans son intégralité.

\smallskip

Une fois que ces paramètres sont connus il est possible de revenir à l'analyse \og patient par patient\fg. En effet,
la formule \eqref{patient} décrit la probabilité $P(D_i \mid \theta_i)$ d'observer les mesures $(d_{i,l})_l$ chez le patient~$i$
connaissant ses paramètres individuels. En utilisant à nouveau la formule de Bayes il vient
\beq \label{Bayespatient}
P(\theta_i \mid D_i)= Z_i^{-1} P(D_i \mid \theta_i)
\eeq
qui fournit la loi de probabilité suivie par le paramètre individuel~$\theta_i$. À~nouveau cette loi peut être approchée
par son mode le plus probable, à savoir le jeu de paramètres $\theta_i$ qui maximise \hbox{$P(\theta_i \mid D_i)$}.

À noter que si il n'y a aucune donnée disponible pour ce patient, alors la formule \eqref{patient} indique
que $\theta_i =\theta^m$: si rien n'est connu sur le patient, il est naturel de supposer qu'il se comporte comme un patient \og moyen\fg.
Au contraire, plus il y a de données sur ce patient, plus le second produit dans \eqref{patient} prend de l'importance.
Lorsque le nombre de données croît, $\theta_i$ se rapproche du jeu de paramètres qui minimise~$J_i$.

\smallskip

Il est maintenant possible de tenter de prévoir l'évolution du patient $i$. Pour commencer il est naturel de supposer que
ses paramètres individuels sont $\theta_i$. Le modèle prévoit que la taille de la tumeur au temps $t$ sera $F(t,\theta_i)$, ce qui
se calcule aisément. Pour ajouter une variabilité sur la prédiction on peut considérer que $\theta_i$ suit la
loi ${\cal N}(\theta^m,\theta^\sigma)$. Les observations $F(t,\theta_i)$ suivent alors la mesure image de cette gaussienne
par $F$. Ceci fournit des intervalles de confiance et des variances de prédiction.

\smallskip

À noter qu'il existe de nombreux tests statistiques et représentations graphiques pour vérifier la validité de cette approche,
tests que nous ne développerons pas ici. Par exemple il est possible de vérifier
que la différence entre prédiction et observation suit bien une loi gaussienne de moyenne nulle.

\subsection{Modélisation des gliomes}

Cette approche de population a été appliquée au gliomes dans l'article \cite{Ribba} sur des données de gliomes
traités par PCV, TMZ et radiothérapie. Certains paramètres du modèle sont relatifs à la croissance du gliome
$(k_{Q_P P}$, $\lambda_P$, $k_{PQ}$ par exemple) alors que d'autres dépendent de la thérapie ($\gamma_P$,
$\gamma_Q$, $\delta_{QP}$). Il est donc possible de traiter toutes les données simultanément, qu'elles viennent
de PCV, de TMZ ou de radiothérapie, en supposant que les paramètres du modèle de croissance sont communs à
toutes ces données, et qu'au contraire les paramètres $\gamma_P$, $\gamma_Q$, $\delta_{QP}$ dépendent de la
thérapie. Cela augmente le nombre de composantes de $\theta$ de $12$ quantités (moyennes et écarts types
pour ces trois paramètres pour le TMZ et la RT), mais augmente significativement le nombre de données patients
utilisables.

L'algorithme SAEM fournit les moyennes et écarts types de tous les paramètres de population.
Il est ensuite possible d'étudier les patients individuellement, en suivant la méthode présentée précédemment.
La~figure \ref{figure7} montre quelques résultats, trois patients pour chaque thérapie.
Les résultats sont qualitativement bons. La courbe noire (paramétrisation du patient) est proche des observations
et décrit bien la dynamique. Les test statistiques sont aussi tous positifs.

\begin{figure}[htb]
\centerline{\includegraphics[width=10cm]{xups22-02_figures/figure7.png}}
\caption{\small Approche populationnelle pour les \hbox{traitements} par PCV, TMZ et radiothérapie (RT). \hbox{Observations}~(cer\-cles),
prédictions individuelles (ligne noire) et intervalles de confiance (zone grisée), prédictions utilisant le paramètre
individuel moyen (en pointillé). }
\label{figure7}
\end{figure}

\smallskip

Il est maintenant temps de s'interroger sur ce qui vient d'être fait.
Nous sommes partis d'un scénario biologique basé sur des hypo\-thè\-ses (existence de deux types de cellules, les proliférantes
et les quiescentes, hypothèses sur le rôle du médicament). Nous avons bâti un modèle d'équations différentielles
ordinaires qui met en jeu ces agents et ces hypothèses. Il n'existe pas une manière unique de mettre en équations
des hypothèses biologiques. Deux modélisateurs différents vont construire deux modèles différents, et d'ailleurs
plusieurs modèles ont été essayés avant de mettre au point le modèle publié.
Ce~modèle a été confronté à des données réelles et a été capable de les reproduire.

Ceci indique que les données cliniques sont cohérentes avec les hypo\-thèses biologiques qui peuvent donc les expliquer.
Il n'y a pas de contradiction. La description \emph{qualitative} donnée par le clinicien peut être rendue \emph{quantitative}
par la modélisation,
ce qui permet d'utiliser les valeurs numériques des données cliniques et non pas seulement un critère Recist de croissance /
stationnaire / décroissance.

\subsection{Optimisation du traitement}

Le modèle mathématique permet de réaliser des expériences, bien sûr impossibles à réaliser sur des vrais
patients. En particulier, il~est possible de changer le rythme d'administration du médicament. Au~lieu de $6$ cycles toutes les
deux mois à une certaine dose, il est possible de simuler la réaction d'un patient à l'administration de $12$ cycles
à demi-dose, ou de $6$ cycles tous les quatre mois, ou n'importe quelle combinaison.

\smallskip

La figure \ref{figure8} montre l'effet du passage d'un cycle de chimiothérapie tous les deux mois en un cycle tous les
neuf mois chez un patient. Dans les deux cas le modèle prévoit une rechute, mais le changement de rythme permet
de gagner presque deux ans avant la remontée de la taille du gliome.
La figure \ref{figure9} montre l'effet combiné du passage de $2$~mois à $6$ mois entre deux cycles et d'un changement
de dose. Le~gain avant rechute est encore plus grand (plus de $7$ ans).

\smallskip

Les prédictions du modèle plaident clairement en faveur d'un changement de rythme d'administration du médicament.
Le passage d'un cycle tous les deux mois à un cycle tous les neuf mois semble permettre d'augmenter le temps avant
recroissance du gliome, et donc potentiellement le temps de survie du patient.

\smallskip

Toutefois, avant de se lancer dans un changement de protocole, il~convient de s'interroger sur ses limitations.
Le modèle a été validé sur une durée de $2$ mois entre deux cycles. Il n'est pas évident qu'il soit
encore valide pour des durées entre deux cycles de $6$ ou $9$ mois. En effet, en particulier sur la figure \ref{figure8},
la taille de la tumeur décroît moins vite lorsque l'on espace les cycles. Comme la taille est plus grande, la probabilité
d'apparition de cellules tumorales résistantes à la chimiothérapie augmente également, et donc le risque
d'une rechute précoce, plus précoce que ce qu'indique le modèle, voire plus précoce que le schéma traditionnel
de $2$ mois. S'il y a des cellules résistantes au début de la chimiothérapie, alors le changement de fréquence des
cycles n'aura probablement que peu d'effet. Il faut donc s'attaquer à la modélisation des mutations, sujet abordé dans
la partie suivante.

\begin{figure}[htb]
\centerline{\includegraphics[width=9cm]{xups22-02_figures/figure8.png}}
\caption{\small Exemple de changement du rythme d'administration du PCV chez un patient. Données expérimentales
(points), rythme classique en vert et nouveau rythme en rose. }
\label{figure8}
\end{figure}

\begin{figure}[htb]
\centerline{\includegraphics[width=9cm]{xups22-02_figures/figure9.png}}
\caption{\small Exemple de changement du rythme et de la dose de PCV }
\label{figure9}
\end{figure}

\section{Apparition de résistances}

Hélas pour les patients de l'étude, la chimiothérapie n'a pas éradiqué la tumeur qui a fini par recroître,
souvent après avoir muté et être passée en \og haut grade\fg (grades $3$ ou $4$). Les cellules tumorales sont alors
devenues résistantes au traitement. La compréhension de ces mutations est fondamentale pour
prévoir l'évolution des tumeurs et optimiser leur traitement.

\begin{figure}[htb]
\centerline{\includegraphics[width=9cm]{xups22-02_figures/figure10.png}}
\caption{\small Trois comportements différents chez des patients traités par TMZ }
\label{figure10}
%\vspace*{-\baselineskip}
\end{figure}

La figure \ref{figure10} montre l'évolution de la taille tumorale chez trois patients différents. Le premier (ID35) voit sa
tumeur décroît longtemps après le traitement, pendant près de trois ans, avant de finalement recroître.
Le second (ID49) a une tumeur qui recroît juste après la fin de la chimiothérapie. Enfin pour le troisième (ID44), la
tumeur recroît même pendant la chimiothérapie. \enlargethispage{\baselineskip}
Ces trois patients illustrent trois scénarios différents et semblent correspondre à l'apparition de résistance
après le traitement, pendant le traitement, voire avant le traitement.


Ceci conduit à la question clinique suivante:

\smallskip

\noindent\emph{Est-ce que des clones résistants sont présents avant
le traitement, ou est-ce le traitement qui induit la résistance?}

\smallskip

\noindent La modélisation mathématique permet d'éclairer un peu cette question (voir par exemple \cite{Ollier} et \cite{these}).

\subsection{Modélisation}

L'idée est de regarder si les données des patients sont mieux expliquées en supposant que la résistance est
présente avant le début du traitement ou en supposant qu'elle apparaît sous l'effet du traitement.

Pour cela nous commençons par construire deux modèles, un modèle de résistance préexistante et un modèle
de résistance acquise. Dans le premier modèle nous supposons que la tumeur est composée de trois types
de cellules: des cellules proliférantes sensibles au traitement, notées $S(t)$, des cellules proliférantes résistantes
au traitement $R(t)$ et des cellules endommagées par la chimiothérapie $D(t)$ (qui meurent progressivement).
Les cellules $S(t)$, et $R(t)$ se multiplient en suivant un modèle logistique. Les cellules $S(t)$ sont attaquées
par la chimiothérapie et deviennent des cellules $D(t)$ qui meurent. On note de plus $C(t)$ la concentration
en TMZ.
Cela conduit au modèle suivant
\begin{align*}
\frac{d S}{dt} &= \lambda S( 1 - \sfrac{P}{P_0}) - \tau_{SD} C(t) S,
\\
\frac{d D}{dt} &= \tau_{SD} C(t) S - \mu_D D,
\\
\frac{d R}{dt} &= \lambda R( 1 - \sfrac{P}{P_0})
\end{align*}
où $P = S + D + R$ est la taille totale de la tumeur. Dans ces équations, $\mu_D$ est la mortalité des cellules
endommagées, et $\tau_{SD}$ mesure l'efficacité de la chimiothérapie. À noter que nous supposons
que $S$ et $R$ ont la même vitesse de croissance, ce qui est naturel car si une des populations allait plus vite que l'autre
elle serait très largement prédominante au début du traitement et l'autre ne serait pas observée.

Dans le second modèle nous ajoutons la possibilité que le TMZ induise la création de cellules résistantes $R_C$
qui pourraient avoir une vitesse de croissance supérieure. Cela conduit à\vspace*{-5pt}
\begin{align*}
\frac{d S}{dt} &= \lambda S( 1 - \sfrac{P}{P_0}) - \tau_{SD} C(t) S,
\\
\frac{d D}{dt} &= \tau_{SD} C(t) S - (\mu_D + \tau_{DR_C}) D,
\\
\frac{d R_C}{dt} &= \lambda ( 1 + \Delta) R_C( 1 - \sfrac{P}{P_0}) + \tau_{DR_C} D,
\\
\frac{d R}{dt} &= \lambda R( 1 - \sfrac{P}{P_0})
\end{align*}
où $P = S + D + R + R_C$. Dans ce système, $\Delta \ge 0$ décrit l'écart de vitesse de croissance entre les cellules
$S$ ou $R$ et les cellules résistantes induites par le TMZ, et $\tau_{DR_C}$ les mutations induites par le TMZ (sur les cellules
endommagées).\enlargethispage{1.5\baselineskip}%

Ces modèles sont alors paramétrés par une approche de population, en ajoutant le fait que chaque patient peut
suivre le premier ou le second modèle. L'algorithme trouve la distribution des paramètres dans la population,
et pour chaque patient, s'il est mieux décrit par le premier ou par le second modèle.

\subsection{Quand apparaît la résistance?}

Le résultat est que pour $51 \%$ des patients, les données cliniques sont mieux expliquées par le second modèle que
par le premier: la chimiothérapie par TMZ provoque elle même l'apparition de résistance chez le patient dans la moitié
des cas. Dans l'autre moitié des cas la résistance au TMZ est déjà présente au moment du traitement.

\begin{figure}[htb]
\centerline{\includegraphics[width=.9\textwidth]{xups22-02_figures/figure11.png}}
\caption{\small Apparition de résistance chez des patients soignés par TMZ: données cliniques (points), taille prédite
par le modèle (courbe noire), $S$ (vert), $D$ (jaune), $R$ (marron), $R_C$ (rouge) }
\label{figure11}
\vspace*{-\baselineskip}
\end{figure}

La figure \ref{figure11} montre l'analyse des données de quelques patients. Les trois figures du haut illustrent
les différentes variables du modèle et les données cliniques chez des patients où le TMZ induit une résistance.
La recroissance de la tumeur est clairement due à la croissance brutale de la courbe rouge (cellules $R_C$) qui deviennent
prédominantes: le~clone résistant induit une aggravation du gliome qui se met à croître beaucoup plus vite
(possible passage en stade $4$).
Les trois figures du bas illustrent des cas où la résistance était déjà présente initialement. La courbe marron
pilote la croissance après disparition des cellules sensibles $S$. La chimiothérapie n'a fait que tuer une petite partie
de la tumeur initiale.

Les patients du premier modèle ont une survie moyenne ($50$ mois) plus faible que ceux du second modèle
($140$ mois), ainsi qu'une survie sans progression inférieure ($22$ mois contre $49$ mois).
Il est aussi possible de corréler l'appartenance à un groupe ou à un autre à l'âge du patient, au profil moléculaire
de la tumeur (sur expression du gène $p53$, mutation IDH ou co-délétion 1p/19q) ainsi qu'à la taille de la tumeur
au début du traitement.

Par exemple, le patient a moins de risque d'appartenir au second modèle en cas de co-délétion 1p/19q (une mutation
classique dans ce type de tumeurs). Ceci est cohérent avec une étude clinique randomisée qui indique que le traitement
par TMZ est meilleur que la radiothérapie en cas de présence de cette délétion, et moins bon dans le cas contraire.

\subsection{Comment optimiser un traitement?}

La biologie sous-jacente aux deux groupes de patients est très différente. Dans un cas les cellules résistantes
sont déjà présentes, dans l'autre elles apparaissent sous l'effet du traitement. Il semble logique que le protocole
thérapeutique optimal soit différent dans les deux cas.

Dans le premier cas, le traitement n'agit que sur une faible partie de la tumeur initiale. Une bonne partie de la tumeur
n'est pas affectée et continue sa progression sous traitement. L'effet du traitement est mineur et il y a peu de place
pour son optimisation.
La~figu\-re~\ref{figure12} montre une tentative a posteriori d'optimisation du traitement d'un patient.
Le protocole $1$ comprend $20$ cycles, tous les $72$ jours à $120mg/m^2/j$ (en cancérologie les doses sont
exprimés en $m^2$ de surface corporelle), le protocole $2$ comprend $10$ cycles tous les $135$ jours à $200 mg/m^2/j$
et le protocole $3$ comprend $15$ cycles tous les $55$ jours à $90 mg/m^2/j$. L'optimisation est a posteriori
car le patient est traité avec le protocole standard, puis ses données sont analysées par le modèle pour en déduire
ses paramètres individuels $\theta_i$, ensuite nous
calculons l'évolution que prévoit le modèle pour ce patient si on change le protocole.
Les trois protocoles conduisent à un gain de survie comparable, de l'ordre de quelques mois.

\begin{figure}[htb]
\centerline{\includegraphics[width=8cm]{xups22-02_figures/figure12.png}}
\caption{\small Tentative d'optimisation de traitement pour un patient du modèle $1$. Les courbes rouge, bleue et rose
correspondent aux protocoles $1$, $2$ et $3$ et la courbe noire au protocole standard. }
\label{figure12}
\end{figure}

La figure \ref{figure13} montre une tentative similaire sur un patient suivant le second modèle.
Les protocoles sont les suivants: $10$ cycles tous les $150$ jours à $200mg/m^2/j$, $15$ cycles tous les $130$ jours
à $200 mg/m^2/j$ et $20$ cycles à $115$ jours à $200 mg/m^2/j$. Le gain en survie prévu par le modèle
est bien plus important, de l'ordre de plusieurs années. Comme la résistance n'est pas présente initialement,
il est possible d'optimiser le traitement pour contrôler $S$ sans créer trop de $R_C$.

\begin{figure}[htb]
\centerline{\includegraphics[width=8cm]{xups22-02_figures/figure13.png}}
\caption{\small Tentative d'optimisation de traitement pour un patient du modèle $2$. Les courbes rouge, bleue et rose
correspondent aux protocoles $1$, $2$ et $3$ et la courbe noire au protocole standard. }
\label{figure13}
\end{figure}

\subsection{Vers une médecine personnalisée?}

Les expériences \og virtuelles\fg du paragraphe précédent indiquent que le traitement par TMZ pourrait être optimisé
et conduire à des gains à la fois en survie et en temps de vie sans progression tumorale qui pourraient être de l'ordre
de quelques mois à quelques années.

L'implémentation d'une telle optimisation sur des patients soulève toutefois de nombreuses questions.
Bien sûr une question éthique tout d'abord: si le modèle est incorrect, le protocole soit disant optimisé pourrait
être pire que le protocole standard, ce qui reviendrait à avoir une action délétère sur le patient.
Ensuite une question fondamentale: l'optimisation est simple mathématiquement si les paramètres individuels
du patient sont connus, mais ce n'est pas le cas en pratique.

\smallskip

Cette seconde question peut être abordée en ayant à nouveau une approche bayesienne.
Pour simplifier, supposons que l'on sache que le patient suive le second modèle.
Soit $\theta$ l'ensemble des paramètres de ce modèle. Au fur et à mesure que le patient passe des IRM,
des nouvelles données cliniques sont disponibles. La formule \eqref{Bayespatient} donne alors la loi des paramètres
individuels $\theta_i$ de ce patient, connaissant ces données.
Au fur et à mesure que les données cliniques s'accumulent, cette loi part
de ${\cal N}(\theta^m,\theta^\sigma)$ pour se rapprocher du paramètre du patient.

\smallskip

La figure \ref{figure14} illustre cette approche, développée dans \cite{these}:
les données des patients jusqu'à trois mois après le début du traitement
son analysées par cette méthode, ce qui fournit une approximation du paramètre individuel $\theta_i$. L'évolution ultérieure
du patient est alors prédite (courbe rouge) et comparée aux données cliniques non incluses dans l'analyse (cercles roses).
La prédiction est assez correcte, en particulier du temps de recroissance. La taille minimale est prédite avec une erreur
inférieure à $1$ cm, ce qui est cliniquement acceptable, dans $90\%$ des cas.
Le paramètre individuel obtenu en utilisant seulement les données cliniques jusqu'à $3$ mois après le début de la
chimiothérapie permet donc de prédire avec une bonne qualité l'évolution du patient et sa rechute.
En utilisant ce paramètre il est alors possible d'évaluer le nombre de cycles de chimiothérapie conduisant à une
rechute la plus tardive possible. La figure \ref{figure15} illustre l'intérêt d'ajouter $5$ cycles chez un patient particulier.
Dans cette étude les simulations du modèle indiquent que
$60\%$ des patients de l'étude auraient bénéficié de cycles de TMZ supplémentaires.
Pour les $40\%$ restants, la résistance au traitement aurait pu être anticipée, et ces patients auraient pu bénéficier
d'un changement de thérapie plus rapidement \cite{these}.

\begin{figure}[htb]
\centerline{\includegraphics[width=9cm]{xups22-02_figures/figure14.png}}
\caption{\small Prédiction de l'évolution sous TMZ à partir des données disponible jusqu'au troisième mois de
traitement. Données intégrées dans la prédiction (bleu), données (rose), prédiction (rouge) }
\label{figure14}
\end{figure}

\begin{figure}[htb]
\centerline{\includegraphics[width=8cm]{xups22-02_figures/figure15.png}}
\caption{\small Bénéfice de l'ajout de $5$ cycles de chimiothérapie chez un patient. Données intégrée
dans l'évaluation des paramètres individuels (rose), prédiction du modèle avec $5$ cycles en plus (bleu) et sans ces cycles
additionnels (bleu puis noir) }
\label{figure15}
\end{figure}

\section{Modèles tri-dimensionnels}

L'IRM ne permet de détecter le gliome que là où sa densité est assez élevée. Si la densité des
cellules tumorales est en dessous d'un certain seuil $\sigma$, elles ne sont pas vues.
Il est possible d'imaginer que la zone visible du gliome diminue, alors que la zone invisible continue de croître.
Le volume total pourrait augmenter alors que par un effet de réorganisation ou un biais de mesure le volume apparent
pourrait diminuer.

\smallskip

Ceci conduit à la question clinique suivante:

\smallskip

\noindent\emph{Et si la décroissance apparente de la tumeur n'était qu'une illusion d'optique?
Et si en fait la tumeur croissait tout le temps, mais que l'IRM donnait l'illusion qu'elle décroissait?}

\smallskip

\noindent Cette question conduit en fait à une question délicate d'équations aux dérivées partielles.

\subsection{Modélisation}

Les modèles ci-dessus sont beaucoup trop compliqués pour pouvoir être étudiés directement. Nous allons donc
présenter un modèle classique beaucoup plus simple.
Soit $u(t,x)$ la densité de cellules tumo\-rales à l'instant $t$ et au point $x$. Les cellules tumorales migrent aléatoirement
et se multiplient, ce qui conduit au modèle dit \og KPP\fg pour Kolmogorov Petrovsky Piskunov ou \og Fisher\fg
\beq \label{Fisher}
\partial_t u - \nu \Delta u = \lambda u (1 - u).
\eeq
L'aspect \og migration\fg est pris en compte par le terme de diffusion, et l'aspect multiplication par un terme logistique.
Le coefficient $\nu > 0$ est un coefficient de diffusion, et $\lambda$ est la vitesse maximale de croissance.
La densité $u$ a été normalisée par la capacité totale du milieu pour varier entre $0$ et $1$.
Pour ce modèle, les cellules sont tout le temps en multiplication, et la masse totale de la tumeur
(hors traitements) est toujours croissante. En effet cette masse vaut
\[
m(t) =\int u(t,x) dx
\]
et en dehors des temps de traitement,
\beq
\dot m(t) = \partial_t \int u(t,x) dx =\lambda \int u (1 -u) dx > 0.
\eeq
Il faut maintenant modéliser le traitement et les IRM.
Pour le traitement, on suppose pour simplifier qu'une fraction $1- \beta$ des cellules est simplement tuée au moment
du traitement, ce qui conduit à
\beq \label{Fisher2}
u(T^+,x) =\beta u(T^-,x)
\eeq
où $T$ est le temps du traitement. Son effet est supposé immédiat: juste après $T$ la densité est multipliée par
$\beta$. Cette hypothèse est assez raisonnable car les molécules de la chimiothérapie disparaissent en quelques
jours alors que la décroissance apparente du gliome se poursuit sur des mois, donc sur une échelle de temps
beaucoup plus grande.

Pour l'IRM, nous faisons l'hypothèse que seules les zones où $u > \sigma$ sont détectées, ce qui conduit à
un volume observé
\beq \label{Fisher3}
V(t)= \int 1_{u(t,x) >\sigma} dx.
\eeq
La question de savoir si la décroissance observée après traitement est un artefact de mesure peut alors se reformuler
ainsi: est il possible que le volume observé $V(t)$ décroisse, alors que la masse tumorale $m(t)$ croît toujours?
En particulier est il possible d'avoir $V'(T^+) < 0$?

\subsection{Étude mathématique}

L'équation \eqref{Fisher} est une équation dite de réaction-diffusion, car elle combine deux phénomènes: une diffusion
(marche aléatoire, diffusion thermique), modélisée par le laplacien, et une réaction (réaction chimique, multiplication,
reproduction,...) modélisée par le terme non linéaire $\lambda u(1-u)$, ou plus généralement $f(u)$, $f$ étant
une fonction quelconque.

Cette équation apparaît très naturellement dans de nombreux domai\-nes scientifiques.
Par exemple en écologie lorsque l'on étudie l'envahissement d'un espace par une nouvelle espèce, il faut étudier
l'interaction entre les mouvements aléatoires des animaux (modélisés par un laplacien) et leur multiplication
(modélisée par un terme logistique). En chimie, la combinaison entre diffusion de réactants et réactions chimiques
conduit à des flammes ou ondes d'explosion, elles aussi modélisées par ce type d'équations. Cette équation
modélise aussi la propagation de fronts d'épidémies (le lecteur intéressé est renvoyé à l'excellent livre de J.D.\,Murray
\cite{Murray}).

Cette équation est aussi très étudiée en mathématiques depuis le grand mathématicien russe Andreï Kolmogorov.
La principale carac\-té\-ristique de ces équations est que le comportement des solutions en temps grand est
complètement connu. Plus précisément, soit $u_0(x)$ une fonction donnée régulière (disons $C^\infty$)
à support compact telle que $0 \le u_0(x) \le 1$ pour tout $x \in \rit^n$. Alors il existe une unique solution
$u(t,x)$ à \eqref{Fisher}. Cette solution est $C^\infty$ en ses deux variables $t \ge 0$ et $x \in \rit^n$.
De plus $0 \le u(t,x) \le 1$ pour tous $t$ et $x$.

Le comportement asymptotique est connu. Tout~d'abord, $u(t,x)$ tend vers $1$ quand $t \to +\infty$,
uniformément sur tout compact. Plus précisément, $u(t,x)\!\sim\!U(\| x \| - ct)$, où $U$ est une fonction explicite et où
$c =\sqrt{2 \nu \lambda}$ (nous ne préciserons pas ici la signification exacte de~$\sim$). Autrement dit, asymptotiquement,
la solution est radiale et, le long d'un rayon, $u$ passe de $0$ (loin de l'origine) à $1$, suivant un profil $U$ qui se déplace
à vitesse constante. Ce profil décrit la transition entre la zone non tumorale et la zone tumorale, ou entre la zone non
envahie par une espèce et la zone envahie.

Revenons à notre question.
Pour calculer le signe de $V'(T^+)$, nous commençons par calculer $\partial_t u(T^+,.)$. Il vient
\[
\begin{aligned}
\partial_t u&(T^+,.) =\nu \Delta u(T^+.) + \lambda u(T^+,.) \bigl( 1 - u(T^+,.) \bigr) \\
&= \nu \beta \Delta u(T^-,.) + \lambda \beta u(T^-,.) \bigl( 1 - \beta u(T^-,.) \bigr) \\
&= \nu \beta \Delta u(T^-,.) + \lambda \beta u(T^-,.) \bigl( 1 - u(T^-,.)\bigr)
+ \lambda (\beta- \beta^2) u(T^-,.)^2 \\
&= \beta \partial_t \Delta u(T^-,.) + \lambda (\beta -\beta^2) \Delta u(T^-,.)  \ge \beta \partial_t u(T^-,.).
\end{aligned}
\]
Autrement dit, si on a $\partial_t u(T^-,x) \ge 0$ avant le traitement, alors, juste après le traitement, on a $\partial_t u(T^+,x) \ge 0$.
Dans ce cas le domai\-ne $\{ x | u(t,x) \sigma \}$ est croissant au sens de l'inclusion en $T^+$ et donc $V'(T^+) \ge 0$.\enlargethispage{\baselineskip}

\smallskip

La question peut alors se reformuler ainsi: est il vrai que toute solu\-tion de \eqref{Fisher} est croissante en temps,
au moins au bout d'un certain temps?

\smallskip

De manière bizarre, si le comportement asymptotique des solutions de~\eqref{Fisher} est bien connu, la question de savoir
si la dérivée temporelle est strictement positive partout au bout d'un certain temps était complètement ouverte.
Il se trouve que cette propriété est vraie, mais sa preuve est très technique \cite{Hamel}.
Le résultat est le suivant: soit $u(t,x)$ une solution de \eqref{Fisher} telle que $0 \le u(0,x) \le 1$
et telle que $u(0,x)$ soit à support compact. Alors il existe $T > 0$ tel que pour tout $t \ge T$, on ait $\partial_t u(t,x) \ge 0$
partout en $x$.

Noter que ce résultat ne peut pas être vrai pour $T = 0$, car si on part de $u_0 = 1_{B(0,1)}$ alors
$u(t,x)$ décroît sur le bord de cette boule à $t = 0$. De façon transitoire, $\partial_t u$ peut être négative,
mais cela ne concerne que le temps initial, donc bien avant que la tumeur ne soit détectée.

\section{Pour aller plus loin}

Sur le sujet de la modélisation du gliome, le lecteur est renvoyé à la thèse de doctorat de P. Mazzocco \cite{these}
pour plus de détails et une bibliographie détaillée. Cette présentation en est largement inspirée.

Les livres de James D.\,Murray sont une excellente introduction à
la modélisation mathématique en biologie et écologie et couvrent une grande variété de sujets, parfaitement
illustrés par des exemples concrets (écologie, morphogenèse...).
Le chapitre $11$ du volume $2$ est en particulier consacré à la modélisation tridimensionnelle de tumeurs cérébrales.

La modélisation mathématique en cancérologie est un sujet en plein développement. En France, de nombreuses
équipes travaillent sur ce sujet, en particulier à l'INRIA (voir le site internet d'INRIA).
La bibliographie qui suit est limitée aux travaux utilisés dans cet article et est tout sauf exhaustive.
Toutefois le lecteur intéressé trouvera facilement des équipes par exemple à Paris, Bordeaux ou Marseille
travaillant sur les chimiothérapies métronomiques, l'étude de la dynamique des métastases,
la croissance de métastase pulmonaires ou l'optimisation de thérapies.\enlargethispage{\baselineskip}

\Subsubsection*{Remerciements}

Les auteurs souhaitent en particulier remercier F.\,Ducray (Centre Léon Bérard) à l'initiative de ces travaux et
P.\,Mazzocco (cette présentation s'inspire directement de sa thèse de doctorat). Ces travaux ont été réalisés
en collaboration avec E.\,Ollier.

\Changelback
\backmatter
\bibliographystyle{jepalpha+eid}
\bibliography{xups22-02}
\end{document}
