# Mathematical modeling and data assimilation for biomedical problems.

Annabelle Collin

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

## To cite this version:

Annabelle Collin. Mathematical modeling and data assimilation for biomedical problems. Mathematics [math]. Université de Bordeaux, 2023. tel-04399487

HAL Id: tel-04399487

https://inria.hal.science/tel-04399487

Submitted on 17 Jan 2024

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_2.png" /></figure>
Thèse présentée en vue d'obtenir le diplôme

Habilitation à Diriger des Recherches

École doctorale de Mathématiques et d'Informatique de l'Université de Bordeaux
Spécialité Mathématiques Appliquées et Calcul Scientifique

par

Annabelle COLLIN

MODÉLISATION MATHÉMATIQUE ET ASSIMILATION DE DONNÉES POUR
DES PROBLÈMES BIOMÉDICAUX.

Soutenue publiquement le 7 décembre 2023, après avis favorable de :

<table><tr><td>Roberto NATALINI</td><td>Directeur de recherche, CNR, Rome, Italie</td><td>Rapporteur</td></tr><tr><td>Mazen SAAD</td><td>Professeur des Universités, École Centrale de Nantes</td><td>Rapporteur</td></tr><tr><td>Marcela SZOPOS</td><td>Professeure des Universités, Université de Paris</td><td>Rapporteuse</td></tr></table>

devant le jury suivant :

<table><tr><td>Fabien CRAUSTE</td><td>Directeur de recherche, CNRS, Paris</td><td>Examinateur</td></tr><tr><td>Astrid DECOENE</td><td>Professeure des Universités, Université de Bordeaux</td><td>Examinatrice</td></tr><tr><td>Florence HUBERT</td><td>Professeure des Universités, Aix-Marseille Université</td><td>Examinatrice</td></tr><tr><td>Magali RIBOT</td><td>Professeure des Universités, Université d'Orléans</td><td>Examinatrice</td></tr><tr><td>Clair POIGNARD</td><td>Directeur de recherche, Inria, Bordeaux</td><td>Garant</td></tr><tr><td>Hugues LOISEAU</td><td>Professeur des Universités, PH, CHU de Bordeaux</td><td>Invité</td></tr><tr><td>Olivier SAUT</td><td>Directeur de recherche, CNRS, Bordeaux</td><td>Invité</td></tr></table>
<figure><img src="image_3.png" /></figure>
REMERCIEMENTS

Mes premiers remerciements vont à Clair Poignard. Merci de de m'avoir guidée et encouragée ces dernières années. Tes connaissances scientifiques et ton recul remarquable sur de nombreux sujets sont une source d'inspiration pour moi. Les travaux présentés dans ce manuscrit te doivent beaucoup. Je remercie également grandement Olivier Saut. En plus de tes nombreuses qualités scientifiques qui rendent la recherche avec toi enrichissante, je te remercie pour ton accueil chaleureux au sein de l'équipe MONC. Avoir un chef d'équipe aussi bienveillant au début de ma carrière a été très précieux. Je tiens aussi à remercier Thierry Colin – même si nous n'avons pas eu le temps de travailler ensemble très longtemps au sein de l'équipe – pour la confiance que tu m'as accordée lors de mon arrivée sur Bordeaux. Un grand merci également à Dominique Chapelle : chaque jour, je réalise l'importance de ce que tu m'as transmis, en particulier ta curiosité et ta rigueur scientifique.

Un grand merci à Roberto Natalini, Mazen Saad et Marcela Szopos qui ont accepté de rapporter ce manuscrit. Je suis très honorée de l'intérêt que vous avez porté à mon travail. Je remercie également très chaleureusement Fabien Crauste, Astrid Decoene, Florence Hubert, Magali Ribot et Hugues Loiseau d'avoir accepté de siéger dans mon jury. Votre présence me touche beaucoup.

Je tiens également à remercier l'ensemble de mes collaboratrices et collaborateurs qu'elles ou ils soient mathématicien·ne·s, biologistes ou clinicien·ne·s. Ces collaborations sont extrêmement précieuses pour mon travail. Je remercie tout particulièrement Philippe Moireau : c'est à chaque fois un vrai plaisir de travailler avec toi, et je suis enchantée par toutes les pistes que nous avons encore à explorer. Un merci spécial aussi à Mélanie Prague : la question que tu m'as posée il y a quelques années m'a fait découvrir tout un nouveau monde, celui des biostatistiques, et je suis convaincue que cette découverte a beaucoup apporté à ma recherche. Et pour terminer, merci à Christèle Etchegaray. Je prends beaucoup de plaisir dans cette nouvelle collaboration et j'espère qu'elle sera l'occasion que tu m'apprennes plein de choses sur les équations différentielles stochastiques.

L'équipe MONC a la chance de compter de nombreuses personnes brillantes, enrichissantes et pleines de bonne humeur. Je remercie Baudouin, Nicolas, Benjamin, Cédrick, Thibaut, Giorgia, Julien, Emma, Luc, Van Linh, Khaoula, Kylian, Eloïse, Guillaume et Florian pour leur contribution à une ambiance de travail agréable. Je souhaite remercier encore plus chaleureusement celles et ceux qui m'ont fait ou me font l'honneur de vivre l'expérience de la thèse avec moi : Sergio, Pedro, Virginie, Simone, Simon et Audrey. Ce manuscrit vous doit beaucoup, et je considère comme un véritable privilège d'avoir été ou d'être témoin de votre évolution sur ce chemin exigeant mais tellement riche qu'est la thèse.

Je tiens à remercier les collègues de l'Institut de Mathématiques de Bordeaux et d'Inria avec qui j'ai la chance d'interagir régulièrement pour leur contribution à la convivialité de notre cadre de travail. Mes remerciements s'adressent en particulier aux membres de l'équipe Calcul Scientifique et Modélisation pour l'émulation scientifique qu'elles et ils génèrent. Une partie de mon intégration et de mes rencontres s'est réalisée à travers le groupe Femmes & Sciences : je remercie Lisl de m'avoir invitée à le rejoindre. Cette aventure a conduit à la création de la mission parité de l'IMB, dont je suis très satisfaite de l'évolution. Mes remerciements s'adressent aussi à mes collègues de l'ENSEIRB-MATMECA, avec une pensée spéciale pour mes deux collègues de bureau, Héloïse et Natalie, ainsi que pour la super équipe DDRS et les habitué·e·s de la pause midi en salle de détente Matmeca/Informatique. Certains de ces collègues sont même devenus des ami·e·s très proches : un merci tout particulier à Juliette, Mélanie, Mathieu et Nicolas.
Je remercie mes ami·e·s d'ailleurs : Saverio, Hélène, Ludivine, Nicolas, Claire, Aurore, Marion, Dorian, Camille, Max, Caco et Micky. La vie aurait tellement moins de saveur sans vous !

Je remercie ma famille et en particulier mes parents. Leur soutien infaillible est une source de ré-
confort inestimable. Merci aussi à ma sœur Camille et à mon frère Antoine : je suis très fière de vous,
et votre soutien m'est également très précieux. Je remercie également ma belle-famille, avec une pensée
particulière pour mes nièces et mon neveu, de m'avoir accueillie au sein de votre famille ainsi que pour
l'intérêt que vous portez à mon travail.

Je dédicace ce travail à Lise. Te regarder grandir depuis ces cinq dernières années est la plus belle expérience de ma vie. Mes derniers remerciements sont pour Charles. Je n'ai (toujours !) pas les mots pour exprimer tout ce que tu apportes à ma vie. Les épreuves que nous avons traversées ces dernières années n'ont fait que renforcer le sentiment qu'avec toi, je suis à ma place dans ce monde. Vivement la suite !
À Lise,
<figure><img src="image_4.png" /></figure>
# Contents

Résumé (en français) 1

Introduction 7

A Modeling aspects 7

A.1 Mathematical models 7

A.2 Model order reduction 8

A.3 Theoretical investigation 8

A.4 Numerical approximation 8

B Data assimilation aspects 8

B.1 Problem setting 9

B.2 General strategy 10

B.3 Sequential approaches 11

C Main contributions 16

C.1 Contributions to cardiac electrophysiology 16

C.2 Contributions to tumor growth 17

C.3 Contributions to electroporation modeling 19

1 Contributions to cardiac electrophysiology 21

1.1 Microscopic bidomain model 22

1.1.1 Definitions and notation 22

1.1.2 General model of ion motion 23

1.1.3 Electroneutrality, transmembrane equations and summary of equations 24

1.1.4 Bidomain model simplifications 25

1.1.5 Ionic models and mechanical effects 26

1.2 Mathematical analysis of the periodic homogenization 26

1.2.1 Analysis of the microscopic bidomain model 27

1.2.2 Homogenization of the bidomain equations 30

1.3 Effects of mechanical deformations on the bidomain model 32

1.3.1 Macroscopic bidomain incorporating the mechanical deformations 32

1.3.2 Numerical illustrations 34

1.4 Contributions to data assimilation aspects 36

1.4.1 Context: a state observer using shape derivatives 37

1.4.2 A new state observer based on the topological gradient 39

1.4.3 Application to real data 39
<table><tr><td>2</td><td>Contributions to tumor growth modeling</td><td>43</td></tr><tr><td>2.1</td><td>Advection type models</td><td>44</td></tr><tr><td>2.1.1</td><td>General presentation</td><td>44</td></tr><tr><td>2.1.2</td><td>Free growth tumor model</td><td>45</td></tr><tr><td>2.2</td><td>Illustration on a large cohort of meningiomas</td><td>47</td></tr><tr><td>2.2.1</td><td>Meningiomas</td><td>47</td></tr><tr><td>2.2.2</td><td>Clinical data</td><td>48</td></tr><tr><td>2.2.3</td><td>Descriptive power</td><td>49</td></tr><tr><td>2.2.4</td><td>Classification</td><td>49</td></tr><tr><td>2.2.5</td><td>Predictive power</td><td>52</td></tr><tr><td>2.3</td><td>Tumor model with time-evolving shapes</td><td>54</td></tr><tr><td>2.3.1</td><td>Tumor model and analysis</td><td>54</td></tr><tr><td>2.3.2</td><td>State observer: design and study</td><td>55</td></tr><tr><td>2.3.3</td><td>And in practice?</td><td>57</td></tr><tr><td>3</td><td>Contributions to electroporation modeling</td><td>63</td></tr><tr><td>3.1</td><td>Membrane dynamics</td><td>64</td></tr><tr><td>3.1.1</td><td>Bioimpedancemetry</td><td>65</td></tr><tr><td>3.1.2</td><td>Membrane modeling</td><td>71</td></tr><tr><td>3.2</td><td>Tissue modeling</td><td>77</td></tr><tr><td>3.2.1</td><td>Comparison and calibration of different electroporation models</td><td>78</td></tr><tr><td>3.2.2</td><td>Floating potential</td><td>84</td></tr><tr><td>3.2.3</td><td>Conclusion</td><td>88</td></tr><tr><td>3.3</td><td>Growth of tumor spheroids subjected to pulsed electric field</td><td>89</td></tr><tr><td>3.3.1</td><td>Biological experiments</td><td>89</td></tr><tr><td>3.3.2</td><td>Modeling of volume evolution</td><td>90</td></tr><tr><td>3.3.3</td><td>Modeling the 3D evolution of spheroid</td><td>92</td></tr><tr><td>3.3.4</td><td>Population-based estimation for PDE systems</td><td>96</td></tr><tr><td>4</td><td>Research project</td><td>103</td></tr><tr><td>4.1</td><td>Tumor growth modeling</td><td>103</td></tr><tr><td>4.2</td><td>Electroporation modeling</td><td>104</td></tr><tr><td>4.3</td><td>Data assimilation aspects</td><td>106</td></tr><tr><td colspan="2">My publications (peer-review journals and proceedings)</td><td>107</td></tr><tr><td colspan="2">Bibliography</td><td>111</td></tr></table>
# Résumé (en français)

Ce document présente un résumé de mes travaux depuis mon arrivée à Bordeaux en 2015 en tant que maîtresse de conférences (Bordeaux INP, ENSEIRB-MATMECA). Côté recherche, je suis membre de l'équipe MONC (Modélisation en ONCologie) du Centre Inria de l'Université de Bordeaux ainsi que de l'équipe Calcul Scientifique et Modélisation de l'Institut de Mathématiques de Bordeaux.

Mes recherches portent sur les mathématiques appliquées au domaine biomédical. Mon objectif est le développement d'outils mathématiques innovants qui ont pour but (1) d'aider les biologistes à comprendre un phénomène et (2) d'aider les cliniciennes et cliniciens à établir un diagnostic ou à améliorer le suivi des patientes et patients.

Le manuscrit – rédigée en anglais à l'exception de ce résumé – est composé d'une longue introduc-
tion et de 4 chapitres. Dans l'introduction, je fais une présentation de la stratégie générale que j'utilise
dans ma recherche. La première étape de mon approche consiste à définir des problèmes biologiques ou
médicaux en consultation avec des expertes et des experts de ces disciplines. Ensuite, si ces problèmes
peuvent être résolus ou du moins compris à travers mes domaines d'expertise – que sont la **modélisation
mathématique avec des équations différentielles ordinaires (EDO)** ou **partielles (EDP)**, leurs **études
théoriques et numériques** et l'**assimilation de données** – une collaboration s'ensuit. Mes collaborations
– dont les principales sont listées à la fin de ce résumé – avec des biologistes et des cliniciens sont très
précieuses pour mon travail. Enfin et dans le but d'adopter une approche la plus complète possible, je
développe des outils mathématiques innovants – principalement basés sur l'étude théorique ou numé-
rique des EDP, l'analyse asymptotique et les observateurs – qui permettent de répondre à la fois à des
**questions théoriques** et à des **applications pratiques**. La figure 1 résume cette approche.

Plus précisément, une première partie de mon travail porte sur des aspects théoriques liés à la mo-
délisation ou à l'assimilation de données et une deuxième partie est très orientée vers les applications,
notamment lorsqu'il s'agit de travailler avec des données réelles. Une partie de mon travail de recherche
s'appuie sur des données médicales qui ont les spécificités (1) de contenir des données difficiles à intégrer
dans des modèles mathématiques et (2) qu'une population d'individus (patients/patientes et expériences
biologiques) avec la même pathologie est souvent disponible. Ces spécificités sont suffisamment impor-
tantes pour avoir des impacts sur ma recherche. En plus de collaborer avec des experts biomédicaux,
ma recherche implique également la collaboration avec des mathématiciennes et mathématiciens, dont
la liste principale est donnée à la fin de ce résumé. Au sein de ces collaborations, l'une de mes contri-
butions principales réside dans mon expertise à intégrer des données réelles provenant de problèmes
biomédicaux dans des systèmes EDO et EDP.

À la suite du chapitre introductif, ma recherche est présentée à travers mes trois plus grands champs d'application : l'électrophysiologie cardiaque, la croissance tumorale et l'électroporation. Malgré cette présentation par chapitre, de nombreuses interactions – mises en avant – dans le manuscrit existent.
<figure><img src="image_5.png" /><figcaption>FIGURE 1 – Présentation générale</figcaption></figure>

L'électrophysiologie cardiaque est l'étude de l'onde électrique qui précède la contraction cardiaque. De nombreuses pathologies cardio-vasculaires sont dues à des troubles électro-physiologiques qui perturbent le rythme cardiaque. L'étude de ce phénomène est donc un véritable enjeu de santé publique. Mes principales contributions dans ce domaine d'application – résumées ci-dessous – sont présentées dans le chapitre 1.

(1) Fournir une analyse mathématique de la procédure d'homogénéisation périodique qui conduit au modèle bidomaine macroscopique. Dans ce travail, nous avons comblé une lacune de la littérature en fournissant une analyse mathématique basée sur la convergence double échelle de la procédure d'homogénéisation qui conduit au modèle bidomaine macroscopique, le modèle de référence en électrophysiologie cardiaque. Notre analyse est exhaustive dans le sens où nous fournissons des résultats d'existence et d'unicité, d'adimensionnement et de convergence dans le même cadre mathématique. Ce travail, publié dans [8], est présenté dans la section 1.2.

(2) Étendre le modèle bidomaine microscopique à un modèle bidomaine microscopique prenant en compte les déformations du cœur. Dans ce travail, le couplage complet entre l'électrophysiologie cardiaque et la mécanique cardiaque a été étudié en remontant aux origines du modèle bidomaine afin d'intégrer les effets directement sur les lois de conservation. Ce travail est présenté dans la section 1.1.

(3) Déterminer un modèle bidomaine macroscopique intégrant les déformations. Grâce à notre compréhension de l'homogénéisation théorique des équations bidomaines classiques (contribution (1)) et une homogénéisation formelle, l'impact des déformations mécaniques sur les équations macroscopiques a été déterminé. Les contributions (2) et (3) sont présentées dans la section 1.3 et ont été publiées dans [11].
(4) **Intégrer les données du front de dépolarisation.** Les travaux précédents nous ont permis d'obtenir une compréhension précise du modèle direct. Afin que le modèle soit réellement descriptif voire prédictif pour un patient donné, il est fondamental d'estimer ou de contrôler les incertitudes – des conditions initiales et/ou des paramètres du modèle par exemple – à l'aide des données médicales disponibles pour ce patient. Dans cette partie – présentée dans la section 1.4 – nous avons travaillé sur des données de front de dépolarisation à l'aide de méthodes séquentielles d'assimilation de données. Une partie des résultats obtenus a été publiée dans [5, 12].

Le chapitre 2 présente mes travaux en croissance tumorale. Le développement de modèles décrivant la croissance tumorale offre la possibilité de suivre l'évolution de la maladie chez un patient spécifique, voire de prédire l'efficacité ou non d'une stratégie thérapeutique. C'est pourquoi la modélisation mathématique de la croissance tumorale suscite l'intérêt de la communauté médicale. Mes contributions principales à ce domaine d'application sont résumées dans les trois points suivants.

(5) **Sélectionner des modèles mécanistes pour la croissance tumorale.** Contrairement à l'électrophysiologie cardiaque, où un modèle de référence est établi, il n'existe pas de consensus clair pour la modélisation de la croissance tumorale. Une explication possible réside dans les multiples inconnues du système et les incertitudes persistantes concernant les lois de conservation et de comportement régissant la croissance tumorale. Les modèles utilisés par l'équipe MONC se basent sur des équations d'advection qui décrivent l'évolution temporelle des densités de cellules saines et tumorales (voir section 2.1).

(6) **Développer une stratégie d'estimation des paramètres du modèle qui soit efficace pour les patients atteints de méningiome.** Les modèles de croissance tumorale doivent être adaptés à chaque patient afin d'obtenir des simulations prédictives. Dans ce contexte, de nombreuses données médicales (volume, texture, forme des lésions ...) peuvent être utilisées pour estimer les paramètres du modèle. Nous avons développé une stratégie innovante qui combine un modèle réduit *OD* et une stratégie d'estimation basée sur les effets mixtes non linéaires. Cette stratégie – illustrée sur une large cohorte de méningiomes (tumeurs bénignes cérébrales) – permet d'obtenir des prédictions fiables du volume et de la forme des tumeurs. Ce travail – présenté dans la section 2.2 – a été valorisé par deux publications [17, 24].

(7) **Adapter notre stratégie d'estimation pour tenir compte de l'évolution des formes tumorales dans le temps.** La stratégie précédemment mentionnée se révèle inefficace lorsqu'il s'agit de prendre en compte l'évolution temporelle des formes tumorales. Dans un autre travail, nous avons donc développé et mathématiquement analysé un observateur de type Luenberger afin de surmonter cette difficulté. Ce travail – présenté dans la section 2.3 – a été publié dans [20].

Le chapitre 3 met en avant mes travaux sur la modélisation de l'électroporation. L'électroporation consiste à appliquer des champs électriques pulsés de forte intensité et de courte durée sur des cellules à l'aide de cathéters. L'objectif est de créer des défauts, appelés pores, dans la membrane cellulaire. L'électroporation réversible permet l'introduction de molécules non perméables dans les cellules vivantes sans les tuer, tandis que l'électroporation irréversible conduit à la mort des cellules de la région ciblée. Mes contributions mathématiques dans ce domaine d'application sont résumées ci-dessous.

(8) **Utiliser la modélisation mathématique afin d'améliorer notre compréhension du processus de perméabilisation de la membrane.** Notre objectif est de mettre en place des approches mathématiques pour étudier la dynamique des membranes. Dans un premier travail présenté dans la section 3.1.1, nous avons développé une stratégie d'étalonnage robuste pour un montage expérimental à 4 électrodes me-surant l'impédance de myotubes cellulaires soumis à des impulsions électriques. Nous avons également conçu et paramétré un modèle de circuit équivalent pour quantifier le comportement à court terme de la
conductance membranaire quelques secondes après l'impulsion. Dans un second travail présenté dans la section 3.1.2, nous avons élaboré un modèle d'électroporation à l'échelle de la membrane qui intègre la teneur en eau de la membrane et le potentiel électrique grâce à un modèle bidomaine microscopique. Une stratégie numérique innovante et efficace nous a ensuite permis d'effectuer des simulations numériques de la formation de pores sur la membrane. Ces travaux ont été publiés dans [25] et [29].

(9) **Comparer différents modèles d'électroporation à l'échelle tissulaire à travers une confrontation à des données expérimentales.** Dans un premier travail, présenté dans la section 3.2.1, nous avons comparé les solutions de divers modèles d'électroporation provenant de la littérature. Nous avons examiné leur paramétrisation en utilisant des mesures bioélectriques obtenues à partir d'expériences sur des foies de lapin réalisées par Sel et al. [SCB+05]. Dans un deuxième travail présenté dans la section 3.2.2, le choix des conditions aux limites des électrodes passives – très conductrices et minces – a été étudié. En effet contrairement aux électrodes actives qui imposent soit le potentiel, soit le courant (conditions de Dirichlet ou de Neumann), la valeur de l'équipotentielle des électrodes inactives est une inconnue des modèles. Cela conduit à une contrainte non locale du flux le long de l'interface conductrice appelé problème du potentiel flottant. Nous avons développé un modèle pour résoudre ce problème dans le contexte de l'électroporation. Ces deux travaux ont été publiés dans [15, 18].

(10) **Étudier les effets de l'électroporation sur la croissance tumorale.** Les sphéroïdes sont des modèles de culture cellulaire 3D capables de représenter l'environnement et l'architecture complexes des tumeurs. En utilisant un modèle mathématique avancé, qui s'inscrit dans la modélisation présentée dans le chapitre 2, et en employant des méthodes d'estimation appropriées, nous avons caractérisé la croissance de sphéroïdes après leur exposition à un champ électrique pulsé. Les résultats de cette étude, présentés dans la section 3.3, ont été publiés dans [21].

Un axe de recherche transverse à ces domaines d'application est le développement d'algorithmes pour l'assimilation de données. Plus précisément, nous avons travaillé sur l'utilisation des données d'une population d'individus pour estimer (entre autres) les paramètres d'un système EDO ou EDP. L'exploitation des connaissances et des données que nous avons sur une population suivant une dynamique similaire permet de compenser les incertitudes inhérentes à un individu spécifique. C'est le principe des modèles mixtes qui se sont révélés plus robustes face aux erreurs de mesure et aux minimums locaux dans de nombreux cas d'application. Cependant, cette méthode d'estimation est trop souvent coûteuse en temps de calcul pour être utilisée en pratique sur des systèmes d'EDP. La stratégie que nous avons développé est basée sur une vraisemblance intégrant l'ensemble de la population qui est maximisée à l'aide d'un filtre de Kalman, voir [23]. Nous avons également proposé une version d'ordre réduit de ce filtre adaptée aux systèmes de grande taille. Pour évaluer les performances de l'algorithme (voir Algorithme 6 dans l'intro-duction), nous l'avons appliqué à un modèle pharmacocinétique. Cette méthode permet de plus d'inclure et d'estimer facilement des erreurs de modélisation. Cette stratégie a été utilisée pour déterminer les associations entre les interventions mises en place par le gouvernement français lors de la pandémie du COVID-19 (confinements, couvre-feux, fermetures d'écoles etc.) et les taux de transmission en se basant uniquement sur les données d'hospitalisation. Le défi était de travailler sans connaissance préalable de la dynamique des taux de transmission. Les résultats obtenus ont été publiés dans [22]. Enfin, dans un troisième travail [28], j'ai proposé une stratégie compatible avec les systèmes d'EDP. Les erreurs provenant des conditions initiales sont contrôlées par un observateur de Luenberger qui corrige l'état, tandis que les paramètres sont estimés à l'aide d'un filtre de Kalman réduit à l'espace paramétrique. L'efficacité de la stratégie a été évaluée sur des données d'électroporation de sphéroïdes tumoraux. Ces travaux sont présentés dans le chapitre 3 de cette thèse.
Le chapitre 4 conclut ce manuscrit en présentant brièvement mes projets de recherche à court et moyen termes. Je commence par aborder mes perspectives de recherche concernant la croissance tumorale, en mettant en évidence deux nouveaux axes de recherche. Le premier axe consiste à approfondir notre compréhension du processus métastatique en étudiant les cellules tumorales circulantes. Le deuxième axe consiste à développer une nouvelle collaboration axée sur l'utilisation du propranolol, un médicament initialement utilisé pour des troubles cardiaques, mais qui montre des résultats prometteurs dans le traitement de certaines tumeurs vasculaires. Ensuite, je décris mes perspectives de recherche dans le domaine de l'électroporation. Dans le cadre du projet ANR MIRE4VTach que je porte, je souhaite m'intéresser à l'électroporation pour l'ablation cardiaque. C'est l'occasion rêvée pour moi de coupler deux de mes champs d'expertise. Enfin, dans une dernière partie, j'aborde brièvement mes perspectives concernant le développement d'algorithmes pour l'assimilation de données. L'idée est de continuer à explorer les stratégies – combinant les aspects populationnels et les méthodes séquentielles – que nous avons développées.

## Principales collaborations

Le travail présenté dans cette HDR est le fruit de collaborations. Mes principaux collaboratrices et collaborateurs sont listés ci-dessous :

**Mathématiciennes/Mathématiciens :**

C. Etchegaray, C. Poignard, O. Saut (Inria Bordeaux, MONC), S. Impériale, P. Moireau (Inria Saclay, M3DISIM), M. Prague (Inria Bordeaux, SISTM)

Doctorants/Doctorantes : T. Kritter (2015-18, 20%), A. Gérard (2015-19, 10%), S. Corridore (2016-20, 50%), P. Jaramillo-Agayo (2019-23, 50%), V. Montalibet (2021-..., 50%)

**Biologistes :**

G. Caluori (électrophysiologiste, IHU Lyric Bordeaux), T. García-Sánchez (Universitat Pompeu Fabra, Barcelona), J. Kolosnjaj, M. Golzio, M.-P. Rols (Université de Toulouse, CNRS, laboratoire IPBS), L. M. Mir (CNRS, Gustave Roussy, Université Paris-Saclay)

**Cliniciens :**

J. Engelhardt, H. Loiseau (neurochirurgiens, CHU Bordeaux), P. Jaïs (cardiologiste, IHU Lyric Bordeaux)

## Principales contributions

Mes publications publiées dans des journaux et des actes de conférence avec comité de relecture sont listées en page 109. Des liens sont fournis pour accéder librement aux publications. Pour faciliter la lecture, les références chiffrées de ce manuscrit correspondent à mes publications et les références lettrées à la littérature. D'autres informations sont disponibles dans mon CV détaillé.
<figure><img src="image_6.png" /></figure>
# Introduction

In this chapter, I first present the general strategy I follow in my research, which is mainly focused on the biomedical field, and then my main contributions. The first step in my approach is to define biological or medical problems in consultation with experts in those disciplines. When these problems can be solved or at least understood through my areas of expertise – mathematical modeling with differential equations (ODE or PDE), theoretical and numerical study of these models, and data assimilation – collaboration ensues. Aiming for the most complete approach possible, I am developing innovative mathematical tools based mainly on the theoretical or numerical study of PDEs, asymptotic analysis, and data assimilation (including observers)– that can address both theoretical and practical aspects. In particular, a first part of my work deals with theoretical aspects of modeling or data assimilation. Another important part of my work is very application oriented, especially when it comes to working with real data. This chapter also gives me the opportunity to introduce several mathematical notions that will be used in the rest of this manuscript. Section A briefly introduces the modeling aspects, and Section B explains how to deal with the many uncertain quantities that must be specified for the models to be truly descriptive or predictive. Finally, Section C contains a summary of my main contributions.

## A Modeling aspects

### A.1 Mathematical models

Most classical mathematical models are based on the application of conservation principles and constitutive relations relevant to the physical, biological, or medical questions under investigation. A conservation law is usually expressed mathematically as a continuity equation, a partial differential equation that establishes a relationship between the quantity of a set and the "transport" of that set. It states that the quantity received at a point or in a volume can change only by the quantity flowing into or out of the volume. A constitutive equation is a relationship between two physical quantities that is specific to a material or substance and approximates the response of that quantity to external stimuli, usually in the form of applied fields or forces. These constitutive equations in the biomedical context can be physiological or phenomenological. A phenomenological model is based on observations of a system rather than a physical theory.

The choice of a model is a delicate trade-off between the best possible representation of reality, its numerical cost, and its identifiability with the available data. This choice must be made as a function of the ultimate goal. In [11], for example, where the goal was to study the effects of mechanical deformation on cardiac electrophysiology, we come back to the derivation of the equations governing ionic motion in the intracellular and extracellular regions in the conservation laws to ensure that we do not neglect any effect that is comparable to the effects we want to study. In contrast, in [17] we consider a very simple
tumor growth model. One could probably consider a model for tumor growth that is closer to biology, but the major advantage of our model is that it depends on a very small number of parameters. In this case, we chose to keep the model simple to ensure that it could be effectively personalized for each patient using the sparse observations available in routine clinical practice.

## A.2 Model order reduction

Part of my work concerns model order reduction. Model order reduction is a technique for reducing the computational complexity of mathematical models. On the one hand, the strategy can be based on a reduced state space dimension such as 0D integration in [17] or radial coordinates based on the spherical symmetry of the tumour in [21]. On the other hand, theoretically validated approximations to the original model can be written thanks to well-designed asymptotic analyses such as in [18] for a floating potential strategy or in [8] for a homogenization method by 2-scale convergence.

## A.3 Theoretical investigation

By proving the existence and uniqueness of the solution of the model under consideration, it is possible to determine whether the problem is well posed. Although this question is important and regularly addressed in my papers, it is not the main topic of my research. Nevertheless, I cite the hard work done in this article [8] regarding the mathematical proof of the existence and uniqueness of the microscopic bidomain model, a very classical model of cardiac electrophysiology.

## A.4 Numerical approximation

In most cases, the solution is not explicit and is then approximated using classical numerical methods such as finite difference methods, finite element methods, or volume element methods for the state variable. I use these methods rather than developing them myself, although a significant amount of my research time is devoted to writing code for numerical simulations. As we will see in the following section, in many situations we want to estimate the parameters of the model using available observations. Most strategies for estimating model parameters require solving the direct model multiple times. To achieve acceptable computational times, we choose numerical schemes very carefully and we rely on parallelization of the code – especially when working with 3D systems and a population of subjects that we want to couple together – with the goal of handling very large systems.

# B Data assimilation aspects

When using ODE or PDE models in practical applications, a major difficulty is dealing with the many uncertain quantities that must be specified for the models to be truly predictive. These quantities include modeling errors, initial conditions, geometry, or parameters of the model. Fortunately, collected data provide additional information that can be used to circumvent the uncertainties associated with the definition of the dynamical system. The development and the application of appropriate estimation strategies occupy a significant portion of my research time. The main difficulties arise from the fact that observations can be sparse in content – for example, when only one of the unknowns is observed – in space – for example, when only a subdomain or part of the boundary is observed – and also in time. In addition, they can be difficult to compare with the solution of the models and are often noisy.
## B.1 Problem setting

In this subsection, I briefly present the formalism of most of the cases I faced.

**Model formulation** We denote by $u$, the system state and by $\theta$ the vector containing all the parameters of the system. Let $\Omega$ be a bounded and smooth domain of $\mathbb{R}^d$, $d = 1, 2$ or $3$, the state dynamics over the time window $[0, T]$ is modeled by an operator $f$ – eventually nonlinear – such that

$$ \begin{cases} \partial_t u(t, x) = f(t, u(t, x), \theta), & [0, T] \times \Omega, \\ u(0, x) = u_0(x), & \Omega. \end{cases} \tag{1} $$

We assume that it is well-posed meaning that for a given vector $\theta$ of size $N_\theta$ and for a given sufficient regular initial solution $u_0 \in \mathcal{X}$, there exists a unique solution $u \in L^2(0, T; \mathcal{X})$, where $\mathcal{X}$ denotes the state space.

**Observations** In practice, we do not observe $u$ directly, but we do have observations – also called measurements – $y \in \mathcal{Z}$. These observations $y$ are related to $u$ by the measurement process. Denoting by $h: \mathcal{X} \to \mathcal{Z}$ the observation operator, we have

$$ y(t, \cdot) = h(u(t, \cdot)) + e(t, \cdot), \tag{2} $$

where $e$ represents the error inherent to the measurement process. To simplify the presentation, we have chosen to consider continuous observations.

**Discretization** For the sake of simplicity, the presentation is done on the discretized version. Given $N_{\Delta t} \in \mathbb{N}^*$ and $\Delta t = \frac{T}{N_{\Delta t}}$, we define the $N_{\Delta t} + 1$ time points $t_n = n \times \Delta t$ for $n \in [0, N_{\Delta t}]$. Then concerning the spatial discretization when PDE is considered, one could use indifferently finite element, finite volume or finite difference methods and we define $N_{\Delta x}$ the size of the discretized state. The lowercase will be used to denote the spatial discretization of the state-dependent variables. The subscript $n$ designs the time discretization at time $t_n$ of the time-dependent variables.

**Uncertainties** We assume that the initial state $U_0$ and the values of the parameters represented by a vector of size $N_\theta$ denoted by $\theta$ are not perfectly known. We quantify the level of uncertainty on $U_0$ and $\theta$ by assuming that they can be modeled as a Gaussian disturbance of $U_0^\diamond$ and $\theta^\diamond$ means (corresponding to the a priori parts) and $P_U$ and $P_\theta$ covariances operators

$$ U_0 = U_0^\diamond + \zeta_U, \quad \theta = \theta^\diamond + \zeta_\theta, $$

where $\zeta_U$ and $\zeta_\theta$ verify $\zeta_U \sim \mathcal{N}(0, P_U)$ and $\zeta_\theta \sim \mathcal{N}(0, P_\theta)$.

**Measurement error** By denoting $W$ the covariance operator of the observation errors, we assume that the errors follow the centered Gaussian law

$$ E_n \sim \mathcal{N}(0, W), \quad \forall n \in [0, N_{\Delta t}]. $$
**Objective** We combine the uncertainties by defining

$$ \zeta = (\zeta_U, \zeta_\theta) = (U_0 - U_0^\diamond, \theta - \theta^\diamond) \text{ and } P = \begin{pmatrix} P_U & \\ & P_\theta \end{pmatrix}. $$

The objective is to find the values of the uncertainties that maximize the probability of the uncertainties under the observations $\max_{\zeta} P(\zeta | (Y_n)_{n \in [0, N_{\Delta t}]} )$. The cases which interest us here correspond to cases for which it is difficult to directly known the uncertainties $(\zeta_U, \zeta_\theta)$ from the observations. This is why it is very classical to use Bayes theorem

$$ \max_{\zeta} P(\zeta | (Y_n)_{n \in [0, N_{\Delta t}]}) = \max_{\zeta} \left[ P((Y_n)_{n \in [0, N_{\Delta t}]} | \zeta) P(\zeta) \right]. $$

Assuming that all disturbances are Gaussian, the strategy relies on the minimization of the following least-square functional

$$ \begin{aligned} \mathcal{J}(\zeta) &= \frac{1}{2}(\zeta, P^{-1}\zeta) + \frac{1}{2}\sum_{n=0}^{N_{\Delta t}} (E_n, W^{-1}E_n), \\ &= \frac{1}{2}(U_0 - U_0^\diamond, P_U^{-1}(U_0 - U_0^\diamond)) + \frac{1}{2}(\theta - \theta^\diamond, P_\theta^{-1}(\theta - \theta^\diamond)) \\ &\quad + \frac{1}{2}\sum_{n=0}^{N_{\Delta t}} (Y_n - H(U_n^{\{U_0^\diamond + \zeta_U, \theta^\diamond + \zeta_\theta\}}), W^{-1}(Y_n - H(U_n^{\{U_0^\diamond + \zeta_U, \theta^\diamond + \zeta_\theta\}}))), \end{aligned} $$

where $(U_n^{\{U_0^\diamond + \zeta_U, \theta^\diamond + \zeta_\theta\}})_{n \in [0, N_{\Delta t}]}$ is the approximated solution of the dynamics (1) with the initial condition $U_0^\diamond + \zeta_U$ and the parameters $\theta^\diamond + \zeta_\theta$.

**Remark 1.** Some remarks are to be made on this section:

* The discretized version is presented for convenience only.

* Other models for measurement error could be considered, such as proportional models. It is straightforward to include other hypotheses in the least squares function.

* In some cases, we may also want to estimate the covariances $P_U, P_\theta$ or/and $W$. They can be easily integrated into the minimization.

* In many cases, the variables do not follow a Gaussian distribution. However, after a nonlinear transformation, they may approximate a Gaussian distribution. This is the case when they are nonnegative values that can be transformed with a log transform, or bounded values that can be transformed with a logit transform.

* It is possible to take into account in practice the fact that the data are sampled in time or in space.

## B.2 General strategy

Regarding the estimation strategy, I have gained experience with assimilating real data from biomedical problems into EDO and EDP systems. More specifically, I use very classical strategies in some cases when there are no real difficulties, and I develop very efficient strategies to deal with large systems (such as PDE systems), especially when a population of subjects is available.
**Identifiability** The minimization of the functional can be considered as an inverse problem. Inverse problems are usually ill-posed, in contrast to the well-posed problems usually encountered in mathematical modeling. There are two types of identifiability problems. We begin with **structural** identifiability. The task of uniqueness is crucial in uncertainty identification, since it is important that the given data be sufficient to uniquely identify the parameter or initial condition that is being sought. This question leads to the proof of cost function injection. We continue with **practical** identifiability. Perturbations in the data (i.e. W too large) or insufficient priors (i.e. P too large) may lead to an insufficiently robust estimation procedure. Identifiability must be investigated – theoretically, if possible – and, of course, practically.

**Classical methods** There are many methods to solve this type of problem numerically, such as a non-exhaustive list: gradient descent based methods, stochastic strategies such as Monte Carlo based strategies, or expectation maximisation algorithms. If we consider small ODE systems with no particular difficulty in identifiability, very classical strategies such as the Levenberg-Marquardt method [Lev44, Mar63] or the standard trust-region method [MS83] work very well, see for example [19]. In [15] we have also used a stochastic algorithm based on an unscented transform [JU04] for a static PDE model.

**Nonlinear mixed-effects model** However, this type of procedure is not sufficiently robust in many situations, such as when measurements are sparse and/or noisy and when there are many local minima. In this case, it is possible to compensate for the uncertainties in a given subject using the knowledge and data we have about a population of subjects of size $N_S$ that follow the same dynamics. More precisely, the nonlinear mixed-effects model consists in pooling all subjects and estimating a global distribution of uncertainties in the population [LW82, Ver97, Lav14, RHM$^{+}$14]. For example, the uncertainties of the parameters of each subject $j \in [1, N_S]$ are divided into two types of uncertainties

$$ \zeta_{\theta}^{j} = \zeta_{\theta}^{\text{pop}} + \zeta_{\theta}^{\text{ind},j}, \qquad (3) $$

where $\zeta_{\theta}^{\text{pop}}$ is common for all the subjects and corresponds to the fixed effect and $(\zeta_{\theta}^{\text{ind},j})_{j \in [1, N_S]}$ are random effects. This strategy was used for example in [17] and [21].

## B.3 Sequential approaches

The computational costs of the strategies presented above increase dramatically as soon as one increases the size of the dynamics (PDE or large ODE systems) and/or the size of the population in the case of mixed effects. An alternative is to use sequential data assimilation methods. In this context, System (1) with uncertainties on the initial condition $\zeta_u$ and on the parameters $\zeta_{\theta}$ is called the target model:

$$ \left\{ \begin{array}{l} \partial_t u(t, x) = f(t, u(t, x), \theta(t)), \quad [0, T] \times \Omega, \\ u(0, x) = u_0^*(x) + \zeta_u(x), \quad \Omega, \\ \theta = \theta^* + \zeta_{\theta}. \end{array} \right. \qquad (4) $$

Sequential methods consist in adding a gain operator (also called a filter) denoted by $g$ that defines a new system – called the observer system –

$$ \left\{ \begin{array}{l} \partial_t \begin{pmatrix} \hat{u}(t, x) \\ \hat{\theta}(t) \end{pmatrix} = \begin{pmatrix} f(t, \hat{u}(t, x), \hat{\theta}(t)) \\ 0 \end{pmatrix} + g(y(t, \cdot), \hat{u}(t, x)), \quad [0, T] \times \Omega, \\ \hat{u}(0, x) = u_0^*(x), \quad \Omega, \\ \hat{\theta}(0) = \theta^*. \end{array} \right. \qquad (5) $$
This filter should be built such that

$$ \|\hat{u}(t, \cdot) - u(t, \cdot)\|_{X} \to_{t \to \infty} 0 \text{ and } \|\theta - \hat{\theta}(t)\|_{\mathbb{R}^{N_{\theta}}} \to_{t \to \infty} 0. $$

Note that the parameters also have dynamics. These methods can be either of the Kalman type [KB61, JU97, Eve09] – which are optimal in the sense of a certain criterion [Ben71] – or either of the Luenberger type [Lue63, Lue71] or either a combination to obtain effective strategies [CMLT09].

**Kalman filters** More precisely, Kalman filters have the advantage that they can be fitted to all data and all models once the minimization is rewritten as a least squares functional. In linear cases, this leads to the well-known Kalman-Bucy filter – also called the Kalman filter – originally presented in [KB61]. This original filter is restricted to a linear – model and observation – operators. Various approximations have been developed for nonlinear cases. The most classical choice of gain for nonlinear configurations is the Extended Kalman Filter (EKF). The idea is to replace the nonlinear mappings by their tangent operators. The Unscented Kalman Filter (UKF) [JU97] is an alternative for state and parameter estimation in nonlinear cases. This newer method was originally introduced to circumvent the major computational drawback of the EKF method, namely, the computation of the tangent operators of the model and observation operators. Indeed, UFK does not require differentiation of the model and the discrepancy with respect to the trajectory, but instead relies on particles surrounding the trajectory to compute the sensitivities of the system to the data and uncertainties.

However, the Kalman filters have a major drawback, namely the time-consuming computation of the full covariance matrix (of size $(N_{\Delta x} + N_{\theta}) \times (N_{\Delta x} + N_{\theta})$). A first strategy is to use a reduced order for the covariance operators. A typical example of reduced order is to restrict the uncertainty space to the parameter space. The main idea behind the reduced order strategy is to consider a SVD decomposition of the covariance matrix

$$ P = LQ^{-1}L^{T}, $$

with $Q$ an invertible matrix of small size $r$ and $L$ an extension operator.

For linear operators, the authors show in [MCLT08] that this decomposition is stable in time and that the equation of the covariance matrix reduces to two systems with admissible computation times. Extensions of these two systems have been developed for nonlinear cases. For example, one can cite the EKF of reduced order (RoEKF) [MCLT08] and the UKF of reduced order (RoUKF) [MC11a]. It should be noted that a specific analysis of the error generated by these nonlinear Kalman filters is required.

**Luenberger filters** An alternative method to the reduced-order strategy is to build a filter that is not based on an optimal criterion. This idea – originally introduced by Luenberger [Lue63, Lue71] and also called nudging in the PDE community [HA76, SS90] – is based on defining a filter that is as simple as possible so that the error between the observed trajectory and the observer system tends to zero. However, these Luenberger observers must be adapted to any dynamics and are not compatible with parameter estimation, since – by definition – the parameters have no dynamics.

**Joint state and parameters estimation** One strategy is to combine the sequential strategies presented below to propose a joint state and parameter estimation. To do this, two different gains – one for the state and one for the parameters – must be defined. A Luenberger filter is considered for the state space, and the small dimension of the parameter space allows the use of an optimal reduced-order filter, typically a RoUKF. More precisely, the mean and covariance operators are computed from the empirical mean and
empirical covariance based on $N_σ$ sample points – the so-called sigma points [JU97]. Different choices of sigma-points can be used in practice, see [JU97, MC11a] for more details. We define for given sets of coefficients $\alpha = (\alpha_j)_{1\le j\le N_σ}$ and particles $(v^{(j)})_{1\le j\le N_σ}$, the empirical mean and covariance by

$$ \mathbb{E}_{\alpha, N_σ}(v_{(\cdot)}) \stackrel{\text{def}}{=} \sum_{j=1}^{N_σ} \alpha^i v^{(j)}, \quad \text{Cov}_{\alpha, N_σ}(v_{(\cdot)}, w_{(\cdot)}) \stackrel{\text{def}}{=} \sum_{j=1}^{N_σ} \alpha_j (v^{(j)} - \mathbb{E}_{\alpha, N_σ}(v_{(\cdot)})) (w^{(j)} - \mathbb{E}_{\alpha, N_σ}(w_{(\cdot)}))^T $$

The full algorithm (spatial discretized version and explicit Euler scheme) reads

### Initialization

Priors on initial conditions: $\hat{U}_0^+$

Priors on parameters: $\hat{\theta}_0^+$

Priors on the covariance matrix of the parameters: $Q_0$

Initial extension matrix: $L_0$

Error covariances: $(W_n)_{n\in[0, N_{\Delta t}]}$

Define sigma points: $e^{(j)}, \forall j \in [1, N_σ]$

#### Loop $\forall n \in [0, N-1]$

##### [1] Sampling

$$ \begin{pmatrix} \hat{U}_n^{(j)+} \\ \hat{\theta}_n^{(j)+} \end{pmatrix} = \begin{pmatrix} \hat{U}_n^+ \\ \hat{\theta}_n^+ \end{pmatrix} + \delta L_n \sqrt{(Q_n^+)^{-1}} e^{(j)}, \forall j \in [1, N_\sigma] $$

##### [2] Prediction (Forward)

$$ \begin{pmatrix} \hat{U}_{n+1}^{(j)-} \\ \hat{\theta}_{n+1}^{(j)-} \end{pmatrix} = \begin{pmatrix} \hat{U}_n^{(j)+} \\ \hat{\theta}_n^{(j)+} \end{pmatrix} + \Delta t \begin{pmatrix} F(t_n, \hat{U}_n^{(j)+}, \hat{\theta}_n^{(j)+}) \\ 0 \end{pmatrix}, \forall j \in [1, N_\sigma] \quad (6a) $$

##### [3] Correction (Analyze)

### State Correction

$$ \begin{pmatrix} \hat{U}_{n+1}^{(j)-} \\ \hat{\theta}_{n+1}^{(j)-} \end{pmatrix} = \begin{pmatrix} \hat{U}_{n+1}^{(j)-} \\ \hat{\theta}_{n+1}^{(j)-} \end{pmatrix} + \Delta t \begin{pmatrix} G(\hat{U}_{n+1}^{(j)-}, Y_n) \\ 0 \end{pmatrix}, \forall j \in [1, N_\sigma] \quad (6b) $$

### Parameter Correction

$$ L_{n+1} = \text{Cov}_{\alpha, N_\sigma}(e_{(\cdot)}, \hat{U}_{n+1}^{(\cdot)+-}) $$

$$ \Lambda_{n+1} = \text{Cov}_{\alpha, N_\sigma}(e_{(\cdot)}, (H(\hat{U}_{n+1}^{(\cdot)+-}) - Y_n)) $$

$$ Q_{n+1} = \text{Cov}_{\alpha, N_\sigma}(e_{(\cdot)}, e_{(\cdot)}) + (\Lambda_{n+1})^T (W_n)^{-1} \Lambda_{n+1} $$

$$ K_{n+1} = L_{n+1} (Q_{n+1})^{-1} \Lambda_{n+1} (W_n)^{-1} $$

$$ \begin{pmatrix} \hat{U}_{n+1}^+ \\ \hat{\theta}_{n+1}^+ \end{pmatrix} = \begin{pmatrix} \mathbb{E}_{\alpha, N_\sigma}(\hat{U}_{n+1}^{(\cdot)+-}) \\ \mathbb{E}_{\alpha, N_\sigma}(\hat{\theta}_{n+1}^{(\cdot)+-}) \end{pmatrix} + K_{n+1} (Y_n - \mathbb{E}_{\alpha, N_\sigma}(\hat{U}_{n+1}^{(\cdot)+-})) \quad (6c) $$
One can notice regarding the algorithm that the time evolution of the parameters is conditioned only by the Kalman filter (6c), which means

$$ \hat{\theta}_{n+1}^{(j)+} = \hat{\theta}_{n+1}^{(j)-} = \hat{\theta}_{n}^{(j)-}, \forall j \in [1, N_{\sigma}], $$

while the state evolves first with the model dynamics (prediction (6a)), then with the state observer (state correction (6b)), and finally by the Kalman filter (parameter correction (6c)). This last step is fundamental to integrate the sensitivity of the parameters to the state and is the key point of the Reduced-order Unscented Kalman Filter (RoUKF) introduced in [MC11b].

In summary, the general strategy for applying joint state and parameter estimation using observers is divided into two steps:

### (1) Construction and mathematical validation of a Luenberger observer assuming the uncertainties are reduced to the initial condition ($\zeta_{\theta} = 0$ in System (4)).

This consists of defining a gain observer for the state – denoted by $\gamma g_u$ – and proving the decrease of the energy of the error model corresponding to the system verified by $\tilde{u} = u - \hat{u}$,

$$ \left\{ \begin{array}{l} \partial_t \tilde{u}(t, x) = f(t, u(t, x), \theta) - f(t, \hat{u}(t, x), \theta) - \gamma g_u(h(u(t, x)), \hat{u}(t, x)), \quad [0, T] \times \Omega, \\ \tilde{u}(0, x) = \zeta_u(x), \quad \Omega. \end{array} \right. \quad (7) $$

The strictly positive constant $\gamma$ – called the parameter gain – balances the impact of the model correction using the data depending on the level of confidence in them. This justifies the *nudging* terminology used in the data assimilation community [LL08]. The gain is therefore typically adjusted by considering the level of noise in data, but also by evaluating the increased stability introduced by the correction. Sometimes an excessively large gain can be counter-productive – see a striking example in [MCLT09] – and in any case should eventually be optimized – see [VLDP03, SB93, ZNLD92].

### (2) Combination of the Luenberger observer with a Kalman filter using for example Algorithm (6).

**Contributions concerning sequential strategies** In my research I have been involved in the development and theoretical investigation (proof of convergence) of Luenberger observers, but also in the application of this type of strategies. More specifically, I have acquired expertise in the development of estimation strategies for front propagation models – such as reaction-diffusion models, eikonal equations, or general hyperbolic/elliptic models – and front position data. Such models are considered in various applications, such as cardiac electrophysiology, wildfire propagation, and tumor growth modeling.

During my PhD, we proposed a Luenberger observer for reaction-diffusion models, validated it mathematically, and illustrated it numerically with synthetic data. Recent works are directly related to this observer. In the context of cardiac electrophysiology, this observer has been applied to real data [12] and a novel correction term based on topological gradients has been integrated to track solutions of complex patterns [5]. In the context of wildfire spread, a version compatible with an eikonal equation was proposed and illustrated with synthetic data in [9] and then with real data in [14]. In this context, a version compatible with a Lagrangian front-tracking model was also developed [13]. Finally, in the context of tumor growth, a Luenberger observer – inspired by our previous work – was developed, mathematically studied (proof of convergence), and numerically illustrated, see [20] for more details.

Another important research result concerning data assimilation is the design of a population-based Kalman filter. Our objective is to propose a strategy to deal with population-based estimation as introduced previously (see Eq. (3)) for large systems as for example PDE. To do so, we first define population-based estimators through maximum likelihood estimation. We then develop an equivalent robust sequen-
<figure><img src="image_7.png" /></figure>

<table><caption>Figure 2 – Top: Estimation of transmission rate with a population-based Kalman filter for the 12 French regions (red) and regression with seasonal weather conditions, the proportion of variants, and government measures (black). Bottom: Estimated associations with transmission rates.</caption><thead><tr><th>Covariate</th><th>Modification of the transmission rate</th></tr></thead><tbody><tr><td><strong>NPIs</strong></td><td></td></tr><tr><td>Lockdown 1 - delay of 7 days</td><td>-78% [-82%; -74%]</td></tr><tr><td>Post lockdown 1 - Phase 1</td><td>-54% [-56%; -52%]</td></tr><tr><td>Post lockdown 1 - Phase 2</td><td>-48% [-50%; -47%]</td></tr><tr><td>Lockdown 2 - delay of 7 days</td><td>-54% [-57%; -49%]</td></tr><tr><td>Lockdown 2 with opened shops</td><td>-51% [-53%; -49%]</td></tr><tr><td>Closing schools</td><td>-7% [-8%; -5%]</td></tr><tr><td>Barrier gestures</td><td>-46% [-48%; -44%]</td></tr><tr><td>Curfew at 6PM</td><td>-30% [-33%; -26%]</td></tr><tr><td>Curfew at 8PM</td><td>-28% [-31%; -25%]</td></tr><tr><td>Bar &amp; restaurant closure (ref. W = 0)</td><td>-10% [-13%; -8%]</td></tr><tr><td>Bar &amp; restaurant closure, summer (vs. ref.)</td><td>-8% [-11%; -4%]</td></tr><tr><td>winter (vs. ref.)</td><td>-11% [-14%; -8%]</td></tr><tr><td><strong>Other factors</strong></td><td></td></tr><tr><td>100% of VoC circulating</td><td>22% [15%; 28%]</td></tr><tr><td>Weather effect during summer (ref. W = 0)</td><td>-22% [-24%; -21%]</td></tr><tr><td>during winter (ref. W = 0)</td><td>10% [9%; 11%]</td></tr></tbody></table>

tial estimator for large populations based on filtering theory that sequentially integrates data. Finally, we limit the computational complexity by defining a reduced-order version of this population-based Kalman filter that clusters subpopulations with common observational backgrounds. The performance of the resulting algorithm is evaluated against classical pharmacokinetics benchmarks. An important advantage of the developed method is that model errors can be easily included and estimated.
In [22], the strategy was used to assess associations between nonpharmaceutical interventions – such as lockdowns, curfews, or school closures – implemented to control the COVID-19 pandemic and transmission rates. The challenge was to work without prior knowledge of the dynamics of transmission rates and to use only hospitalization data, which are the only accurate data available. Figure 2 illustrates the results, which were used by the French government to adjust the measures and also cited in French newspapers, see for example [Ros21, Hec21, Don21].

However, our first goal in developing a population-based Kalman filter was to overcome the challenge of a nonlinear mixed-effects strategy for large systems such as PDE systems. Indeed, the nonlinear mixed-effects strategy is not applicable in practice for PDE systems due to its high computational cost. For this reason, we propose in [28] a population-based estimation for PDE systems. The errors arising from the initial conditions are controlled by a Luenberger observer that corrects for the state, and the parameters are estimated using a reduced-order population-based Kalman filter restricted to the parameter space. This sophisticated method for solving this difficult problem is based on two important existing methodological strategies: (1) the population-based Kalman filter presented in [23], and (2) the joint state-parameter estimation proposed in [CMLT09]. The performance of the resulting algorithm is evaluated on tumor spheroid electroporation (synthetic and real data), see Section 3.3.4 for details.

# C Main contributions

This manuscript presents my research based on my three main areas of interest: Cardiac Electrophysiology, Tumor Growth, and Electroporation Modeling. Despite this chapter-by-chapter presentation, there are many interactions – highlighted in the manuscript – between these areas. However, before going into detail in the following chapters, I provide an overview of the major contributions in this section.

## C.1 Contributions to cardiac electrophysiology

Cardiac electrophysiology is devoted to the study of the electrical wave that precedes contraction of the heart muscle. In Chapter 1, I present my significant contributions in this particular field. The following points serve as a concise summary of these contributions:

(1.1) **Provide a mathematical analysis of the periodic homogenization process leading to the macroscopic bidomain model.** In this work, we fill a gap in the literature by providing a mathematical analysis of the macroscopic bidomain model based on 2-scale convergence. Our analysis is exhaustive in the sense that we provide existence and uniqueness, nondimensionalization, and convergence results within the same mathematical framework. This work, published in [8], is presented in Section 1.2.

(1.2) **Extend the microscopic bidomain model to a microscopic bidomain model that accounts for cardiac muscle deformations.** Our goal is to study the effects of cardiac deformations on cardiac electrophysiology. We go back to the origins of the bidomain model to integrate the effects directly into the conservation laws at the microscopic scale. This contribution is presented in Section 1.1.

(1.3) **Determine a macroscopic bidomain model that integrates the cardiac deformations.** Thanks to our understanding of the theoretical homogenization of classical bidomain equations, we perform a formal homogenization to determine the effects of mechanical deformations on the macroscopic equations. This work is presented in Section 1.3 and has been published (with Contribution 1.2) in [11].
(1.4) **Assimilate depolarization front data.** The previous works have given us an accurate understanding of direct models of cardiac electrophysiology. For the model to be truly descriptive or even predictive for a given patient, we need to control for uncertainties arising from the initial conditions or parameters of the model using the available medical data of this patient. In this part presented in Section 1.4, we worked with depolarization front data and developed an efficient Luenberger observer. This work has been published in two proceedings [5, 12].

LEFT ANTERIOR OBLIQUE VIEW OF LEFT ATRIUM

<figure><img src="image_8.png" /><figcaption>Figure 3 – Depolarization maps of a patient left atrium. Left: Data extracted with Carto® system. Right: Completion of the map using a Luenberger observer integrating depolarization fronts.</figcaption></figure>

**Illustrative example** Our hard work on the bidomain model (Contributions 1.1, 1.2, and 1.3) has allowed us to gain a very good understanding of the direct model and, in particular, whether or not we need to consider mechanical deformations depending on the pathology of interest. A direct application of the observer we developed (Contribution 1.4) is its ability to smooth and complete activation maps. Figure 3 illustrates the potential of our approach on the left atrium of a patient. The left side of the figure shows the measured map of activation times corresponding to the times at which the depolarizing front reaches the area. This map is not smooth and is spatially incomplete due to the difficulty of the recording procedure. Our approach – a combination of a well-chosen PDE system and an efficient data assimilation strategy – allows us to smooth and complete the activation maps, as shown on the right side of Figure 3. A complete and realistic map is very useful for clinicians to identify pathological areas that need to be ablated.

## C.2 Contributions to tumor growth

Chapter 2 presents my work on tumor growth whose the main contributions are:

(2.1) **Select mechanistic models for tumor growth.** Unlike cardiac electrophysiology, where there is a gold standard, there is no such consensus for modeling tumor growth. One explanation is that the unknowns of the system and the conservation laws are still open questions because of the many uncertainties about how tumor growth works. The models used by the MONC team are based on **advection equations** that model the time evolution of the densities of healthy and tumor cells.
(2.2) Develop an efficient estimation strategy to fit a tumor growth model to patients with meningioma. Tumor growth models need to be fitted to each patient to obtain predictive simulations. In this context, many medical data (volume, texture, shape of lesions, etc.) can be used to determine the parameters of the model. We developed an innovative strategy that combines a *OD scale model* and an *estimation strategy* based on nonlinear mixed effects. This strategy has been illustrated on a large cohort of meningiomas (benign brain tumors) providing reliable predictions of their volume and shape. This work – presented in Section 2.2 was highlighted by two publications [17, 24].

(2.3) Adapt this strategy to tumors with time-evolving shapes. The strategy described above is not effective when the shape of the tumor changes over time. In this work, a *Luenberger-type observer* has been *developed* and *mathematically analyzed* to solve this problem. This work – which is presented in Section 2.3 – has been published in [20].

DATA - 4 MRI EXAMINATIONS DURING ONE YEAR

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

RESULTS - USING THE 3 FIRST MRI EXAMINATIONS

<figure><img src="image_10.png" /><figcaption>Figure 4 – Top: MRI examinations of a brain metastasis of a patient (HM Hospitales, Spain) at 4 different time points (given in days) and the 3D segmented tumors. Bottom: Volume description and prediction (left), growth rate estimation (middle) and shape prediction.</figcaption></figure>
**Illustrative example** The accuracy of the model considered, despite its simplicity (contribution 2.1) in combination with a 0D scale model (contribution 2.2) and a Luenberger observer (contribution 2.3), allows us to describe and predict the volume and shape of a tumor. Our strategy is based on our ability to build a spatially distributed growth rate. Figure 4 illustrates the potential of our strategy on a patient brain metastasis. At the top, one can see the data and the segmented 3D tumors. At the bottom, one can see the volume evolution over time with a prediction of the last time point, the estimated spatial distribution of the growth rate, and the prediction of the shape of the tumor.

### C.3 Contributions to electroporation modeling

Chapter 3 presents my work on modeling the electroporation phenomenon. Electroporation is the application of pulsed electric fields of high intensity and short duration to cells using catheters. The goal is to create defects – called pores – in the cell membrane. Reversible electroporation allows the introduction of non-permeable molecules into living cells without directly killing them, while irreversible electroporation leads to the death of cells in the target area. My mathematical contributions to this application area are as follows:

(3.1) **Use mathematical modeling to better understand how the membrane becomes permeable.**
The goal is to develop mathematical approaches to track membrane dynamics. In an initial work presented in Section 3.1.1, we develop a robust calibration strategy for a 4-electrode experimental setup to measure the impedance of cellular myotubes exposed to electrical pulses, and we design and parameterize an equivalent circuit model. This allows us to quantify the short-term behavior of membrane conductance a few seconds after the pulse. In a second work presented in Section 3.1.2, we develop a PDE model of membrane electroporation that combines membrane water content and electric potential through a microscopic bidomain model. An innovative and efficient numerical strategy allows us to perform numerical simulations of pore formation on the membrane. These two papers have been published in [25] and [29].

(3.2) **Compare and improve tissue-scale electroporation models.** In a first work presented in Section 3.2.1, we compare solutions of different electroporation models from the literature and investigate their parameterizations using bioelectrical measurements on rabbit livers performed by Sel et al. [SCB+05]. In a second work presented in section 3.2.2, the choice of boundary conditions for passive electrodes – very conductive and thin – is investigated. In contrast to the active electrodes, which specify either the potential or the current (Dirichlet or Neumann conditions), the value of the equipotential of the inactive electrodes is an unknown factor of the PDE, regardless of the model considered. This leads to a nonlocal flux constraint along the conductive interface, the so-called floating potential problem. We develop a model to solve this problem in the context of electroporation. These two works have been published in [15, 18].

(3.3) **Study the effects of electroporation on tumor growth.** The aim of this work is to study the specificities of the growth of cancer cells in spheroids exposed to a pulsed electric field. Thanks to an advanced mathematical model (within the framework of the modeling presented in Chapter 2) and relevant estimation strategies, our modeling allows characterizing the growth of multicellular spheroids after permeabilization by an electric field pulse. This work, presented in Section 3.3 was published in [21].
SPHEROID ELECTROPORATION

<figure><img src="image_11.png" /><figcaption>Figure 5 – Left: Fits of estimated volumes to measured volumes. Right: PDE simulation of the spheroids. Each column corresponds to a different case: control, electroporation field of 500V.cm⁻¹, electroporation field of 1000V.cm⁻¹ (group B), electroporation field of 1000V.cm⁻¹ (group A).</figcaption></figure>

**Illustrative example** Our acquired knowledge of electroporation modeling (contributions 3.1 and 3.2) combined with a realistic tumor growth model (contribution 2.1) and relevant estimation strategies allowed us to characterize the growth of 3D *in vitro* tumors – called multicellular spheroids – after electroporation. The use of a PDE system allowed us to determine the percentage of cells destroyed and the percentage of cells whose functioning was altered by the effect of electroporation for different values of the electrical field. Figure 5 summarizes the results. It can be seen that for an electrical field close to 1000 V.cm⁻¹ (see group A), the increase of the electrical field seems to boost the resumption of spheroid growth. This is a recurrent problem that also occurs in clinics, where clinicians observe either efficient treatment with reduced tumor growth or tumor relapse [SNdJ+14]. The advantage of our model is that we can determine the trend of tumor growth shortly after electroporation. So, in practice, this model can be used to determine which patients have been treated efficiently and which patients are on the path to relapse. In this last case, this would alert the physician that the treatment was not efficient and that another electroporation cycle is needed to remove the tumor again.
# Chapter 1

# Contributions to cardiac electrophysiology

*The innovative work presented in this chapter has been published in two articles [8, 11] and two proceedings [5, 12].*

## Introduction

The purpose of this chapter is to present my mathematical contributions to cardiac electrophysiology since the end of my PhD. Cardiac electrophysiology is the study of the electrical wave which precedes the cardiac contraction. Modeling of cardiac electrophysiology is a very mature field of research, and the so-called bidomain model [Tun78] is now considered the reference model, while some approximations to it are often considered, in particular the so-called monodomain model [Sac04, CFPT05], and also some specific forms of eikonal equations [Kee91, KS04, KCM$^{+10}$]. These macroscopic bidomain equations are the result of a complex homogenization process starting from a microscopic level associated with the cellular scale, where intracellular and extracellular domains – separated by the cell membrane – can be defined geometrically [NK93, CFS02, RC11].

My mathematical contributions aim to fill in some important gaps with respect to the very classical bidomain model, which can be summarized in the following points:

(1) **Provide a complete mathematical analysis of the periodic homogenization procedure that leads to the macroscopic bidomain model.** The aim of this work is to fill a gap in the literature by providing a complete mathematical analysis based on 2-scale convergence theory of the homogenization procedure that leads to the macroscopic bidomain model. Our analysis is exhaustive in the sense that we provide existence and uniqueness results, nondimensionalization of the equations and 2-scale convergence results – in particular, for the non-linear terms supported on the membrane surface – in the same mathematical framework. The results have been published in [8]. In this work, the cardiac mechanical deformations are neglected.

(2) **Extend the microscopic bidomain model by a microscopic bidomain model that considers deformations of the heart muscle.** Our objective is to examine the complete coupling between cardiac electrophysiology and mechanics, since this has not yet been fully elucidated. More specifically, the effect of mechanics on ionic models – i.e. ordinary differential equations designed to describe the local exchange of ions between different cellular compartments, including through cell membranes – has received considerable attention, see [KNP08, NP04], but in contrast, the effects on the conservation laws (of charges and ions) governing the spatial variations of all quantities of interest have not been
so actively studied before our work. We refer to [CFPS16] for a detailed overview of the latter aspect, see also [CFPS17, CFPS18, FPP$^+$16, QLRRB17]. To do this, it is essential to go back to the origins of the bidomain model to finally clarify the effects of mechanical deformations on it. This work allows to determine a bidomain model at the microscopic level i.e. at the cellular level including the effects of the mechanical deformations.

(3) Determine a macroscopic bidomain model including the effects of the mechanical deformations. Thanks to our understanding of the theoretical homogenization of the classical bidomain equations, we perform a formal homogenization to determine how the effects of mechanical deformations modify the macroscopic equations. The two last points have been published in [11].

(4) Integrate data corresponding to the sparse location of the depolarization front into the macroscopic bidomain model. An accurate understanding of the direct model has been achieved thanks to the previous works. For the model to be truly descriptive or even predictive for a given patient, the uncertainties arising from, among other things, the initial conditions or the parameters of the model must be estimated or controlled using the medical data available for that patient. In this work, we focus on data corresponding to depolarization fronts and we manage them with sequential data assimilation methods. The results have been published in [5, 12].

Since the second contribution allows us to explain the origin of the microscopic bidomain model – even without including the deformations of the myocardium –, we start with this point, see Section 1.1. In Section 1.2, the mathematical analysis of the periodic homogenization procedure leading to the macroscopic bidomain model – corresponding to the first contribution – is presented. In this section, the mechanical deformations are neglected, but it allows to present the homogenization procedure leading to the macroscopic bidomain model, which allows a better illustration of the formal homogenization performed to obtain the macroscopic bidomain model including the heart deformations presented in Section 1.3 corresponding to the third contribution. Some numerical illustrations are included. These simulations show, first, that we are now able to obtain very realistic simulations of cardiac electrophysiology, and second, they allow quantification of the effects of mechanical deformations on cardiac electrophysiology. Based on this very good understanding of the direct model, Section 1.4 – corresponding to the fourth contribution – presents our strategy for dealing with the uncertain quantities that must be preserved when performing model simulations for a given patient.

## 1.1 Microscopic bidomain model

The purpose of this section is to recall the microscopic model – governing ionic motion and electric potentials – underlying the bidomain model and to derive extensions of these equations to account for mechanical deformations. This section is a summary of Section 2 of our article [11]. The novelty arises from the consideration of mechanical deformations.

### 1.1.1 Definitions and notation

The geometric domain occupied by the cardiac tissue at any time $t$ is denoted by $\Omega(t)$, with boundary $\partial\Omega(t)$. We use a Lagrangian description, for which we introduce a reference configuration $\Omega_0$ – not necessarily corresponding to $\Omega(0)$ – and the one-to-one deformation mapping $\phi$ from the reference configuration to the current configuration, giving the position of each material point over time, $\phi: \bar{\Omega}_0 \times [0, T] \rightarrow \mathbb{R}^3$ such that $x(t) = \phi(\xi, t)$. We define the displacement $y$ by $y(\xi, t) = x - \xi = \phi(\xi, t) - \xi$, and the deformation
<figure><img src="image_12.png" /><figcaption>Figure 1.1 – Intracellular and extracellular regions of the heart.</figcaption></figure>

gradient by $\underline{F}(\xi, t) = \vec{\nabla}_{\xi} \phi = \underline{1} + \vec{\nabla}_{\xi} y$. We will assume that $J = \det \underline{F}$ is strictly positive for any $\xi \in \bar{\Omega}_0$, meaning that the deformation mapping does not degenerate anywhere and that $\underline{F}$ is invertible.

In fact, this geometric domain $\Omega_0$ is decomposed into two parts, i.e. $\Omega_0^{(i)}$ the intracellular region and $\Omega_0^{(e)}$ the extracellular region, see Figure 1.1, with the associated regions $\Omega_0^{(i)}(t)$ and $\Omega_0^{(e)}(t)$ in the deformed configuration. Geometrically, $\Omega_0^{(i)}$ and $\Omega_0^{(e)}$ are two connected domains verifying $\bar{\Omega}_0 = \bar{\Omega}_0^{(e)} \cup \bar{\Omega}_0^{(i)}$ and $\Omega_0^{(e)} \cap \Omega_0^{(i)} = \emptyset$. We suppose that the membrane separating the two regions $\Gamma_0 = \partial\Omega_0^{(e)} \cap \partial\Omega_0^{(i)}$ is regular, and we define $\underline{n}_0^{(i)}$ and $\underline{n}_0^{(e)}$ as the unit outward-pointing normal vectors associated with $\Omega_0^{(e)}$ and $\Omega_0^{(i)}$, respectively, with of course $\underline{n}_0^{(i)} = -\underline{n}_0^{(e)}$. Accordingly, $\underline{n}^{(e)}(t)$ and $\underline{n}^{(i)}(t)$ denote the unit normal vectors on the deformed membrane $\Gamma(t)$. In what follows, we will often drop the explicit time dependence in the notation pertaining to the deformed configuration, as in $\Omega_0^{(e)}$ and $\Omega_0^{(i)}$.

### 1.1.2 General model of ion motion

In this section, we give an overview of the derivation of the equations governing ionic motion in the intracellular and extracellular regions, namely, the so-called Nernst-Planck-Poisson equations, see e.g. [Kir10]. We consider ion motion in an isothermal medium with negligible solvent flows. We suppose that we have $N^{(i)}$ and $N^{(e)}$ different species of mobile ions in the intra- and extracellular regions, respectively. The individual conservation gives,

$$ \partial_t c_k^{(\alpha)}|_{\underline{x}} + \vec{\nabla}_{\underline{x}} \cdot j_k^{(\alpha)} = 0, \quad k = 1, \dots, N^{(\alpha)}, \quad \alpha = i, e, \qquad (1.1) $$

where $c_k^{(\alpha)}, j_k^{(\alpha)} = c_k^{(\alpha)} \underline{v}_k^{(\alpha)}$ are the concentration and the flux of the k-th ionic species, respectively, in the given region, and with $\underline{v}_k^{(\alpha)}$ the corresponding velocity. Denoting by $\underline{v}_{k,r}^{(\alpha)}$ the relative velocity of the ions with respect to the substrate, we will use the decompositions

$$ \underline{v}_k^{(\alpha)} = \underline{v}_s + \underline{v}_{k,r}^{(\alpha)}, \quad \underline{j}_k^{(\alpha)} = \underline{j}_{k,s}^{(\alpha)} + \underline{j}_{k,r}^{(\alpha)} = c_k^{(\alpha)}(\underline{v}_s + \underline{v}_{k,r}^{(\alpha)}). $$

The total charge density $\rho^{(\alpha)} = \rho_b^{(\alpha)} + F \sum_k z_k c_k^{(\alpha)}$ – with $\rho_b^{(\alpha)}$ the background charge density, F the
Faraday constant, and $z_k$ the signed number of elementary charges carried by this ion – is related to the electric field $\underline{E}^{(\alpha)}$ by the Gauss law of electrostatics

$$ \vec{\nabla}_{\vec{x}} \cdot \underline{E}^{(\alpha)} = \frac{\rho^{(\alpha)}}{\varepsilon_0 \varepsilon_s} = \frac{1}{\varepsilon_0 \varepsilon_s} \left( \rho_b^{(\alpha)} + \mathrm{F} \sum_{k=1}^{N^{(\alpha)}} z_k c_k^{(\alpha)} \right), \qquad (1.2) $$

with $\varepsilon_s$ the relative permittivity of the medium, and $\varepsilon_0$ the permittivity of free space. In order to determine the flux $J_k^{(\alpha)} = c_k^{(\alpha)} \underline{v}_k^{(\alpha)}$ of the species k, we apply Newton's second law to an ion of this species. The forces undergone by this ion are:

* the friction force $\underline{F}_k^{\text{frict}} = -\frac{\underline{v}_k^{(\alpha)}}{m_k}$, where $m_k$ denotes the so-called mobility of the ion;

* the electric force $\underline{F}_k^{\text{elec}} = z_k q \underline{E}^{(\alpha)}$, with $q$ the charge of a proton;

* the diffusion force (Einstein's diffusion law) $\underline{F}_k^{\text{dif}} = -\frac{k_B T}{c_k^{(\alpha)}} \vec{\nabla}_{\vec{x}} c_k^{(\alpha)}$, where $k_B$ denotes the Boltzmann constant and $T$ the absolute temperature.

Neglecting the inertia and using the equilibrium equation $\underline{F}_k^{\text{dif}} + \underline{F}_k^{\text{frict}} + \underline{F}_k^{\text{elec}} = 0$, we infer the Nernst-Planck equation for the flux $J_k^{(\alpha)} = c_k^{(\alpha)} \underline{v}_k^{(\alpha)}$ of the species k,

$$ J_k^{(\alpha)} = -D_k \left( \vec{\nabla}_{\vec{x}} c_k^{(\alpha)} - \frac{q z_k c_k^{(\alpha)}}{k_B T} \underline{E}^{(\alpha)} \right) + c_k^{(\alpha)} \underline{v}_s, \qquad (1.3) $$

with $D_k = m_k k_B T$ the diffusion coefficient of the species k. The electric field is related to an electric potential $u$ by $\underline{E}^{(\alpha)} = -\vec{\nabla}_{\vec{x}} u^{(\alpha)}$. Using (1.1), (1.2) and (1.3), we finally obtain

Charge conservation : $\partial_t c_k^{(\alpha)}|_{\underline{x}} = -\vec{\nabla}_{\vec{x}} \cdot J_k^{(\alpha)}$, $\forall k \in [1; N^{(\alpha)}]$,

Drift-diffusion flux : $J_k^{(\alpha)} = -D_k \left( \vec{\nabla}_{\vec{x}} c_k^{(\alpha)} + \frac{q z_k c_k^{(\alpha)}}{k_B T} \vec{\nabla}_{\vec{x}} u^{(\alpha)} \right) + c_k^{(\alpha)} \underline{v}_s$, $\forall k \in [1; N^{(\alpha)}]$, \qquad (1.4)

Gauss law : $\Delta_{\underline{x}} u^{(\alpha)} = -\frac{1}{\varepsilon_0 \varepsilon_s} \left( \rho_b^{(\alpha)} + \mathrm{F} \sum_{k=1}^{N^{(\alpha)}} z_k c_k^{(\alpha)} \right)$.

This system is an extension of the Nernst-Planck-Poisson system, with the additional term $c_k^{(\alpha)} \underline{v}_s$ results from the substrate motion.

### 1.1.3 Electroneutrality, transmembrane equations and summary of equations

By examining the orders of magnitude of the physical parameters in the above Gauss law, compared to the time and space scales in consideration, it can be argued – and justified by a mathematical asymptotic analysis, indeed [Mor09, Ric09] – that electroneutrality must prevail away from the membranes, i.e.

$$ \rho_b^{(\alpha)} + \mathrm{F} \sum_{k=1}^{N^{(\alpha)}} z_k c_k^{(\alpha)} = 0, $$
everywhere in the interior of the intra- and extracellular regions. In order to obtain an equation for the electric potential, we need to consider the total electric current and apply to it the electroneutrality assumption. This gives the following equations, see [11] for more details,

$$ \underline{j}_r^{(\alpha)} = -\left(\sigma^{(\alpha)} \vec{\nabla}_{\vec{x}} u^{(\alpha)} + \mathrm{F} \sum_{k=1}^{N^{(\alpha)}} z_k D_k \vec{\nabla}_{\vec{x}} c_k^{(\alpha)}\right) \quad \text{and} \quad \vec{\nabla}_{\vec{x}} \cdot \underline{j}_r^{(\alpha)} = 0. $$

To close the system, the drift-diffusion flux at the membrane interface is modeled: the ions that hit the membrane can be blocked on the surface on each side, or cross the membrane through so-called ionic channels. We denote by $\mu_k^{(\alpha)}$ the molar density per unit surface for species k on either side, and by $I_k^{\text{ion}}$ the molar quantity crossing the membrane from $\Omega^{(i)}$ to $\Omega^{(e)}$, per unit time and per unit surface of the reference configuration. We also introduce the total quantities $I^{\text{ion}} = \mathrm{F} \sum_{k=1}^{N} z_k I_k^{\text{ion}}$ and $\mu^{(\alpha)} = \mathrm{F} \sum_{k=1}^{N} z_k \mu_k^{(\alpha)}$. By considering that conservation dictates that the variation of molar quantity attached to the membrane equals the quantity crossing the membrane from $\Omega^{(e)}$ to $\Omega^{(i)}$ plus the quantity associated with ions coming from $\Omega^{(i)}$ and hitting the membrane, the full system of unknowns $c_k^{(\alpha)}$ and $u^{(\alpha)}$ reads

$$ \left\{ \begin{array}{ll} \partial_t c_k^{(\alpha)}|_{\vec{x}} = -\vec{\nabla}_{\vec{x}} \cdot \underline{j}_k^{(\alpha)}, & \text{in } \Omega^{(\alpha)} \\ \underline{j}_k^{(\alpha)} = -D_k \left(\vec{\nabla}_{\vec{x}} c_k^{(\alpha)} + \frac{qz_k c_k^{(\alpha)}}{k_B T} \vec{\nabla}_{\vec{x}} u^{(\alpha)}\right) + c_k^{(\alpha)} \underline{v}_s, & \text{in } \Omega^{(\alpha)} \\ s_{\Gamma} \underline{j}_{k,r}^{(\alpha)} \cdot \underline{n}^{(\alpha)} = \partial_t \mu_k^{(\alpha)}|_{\xi} \pm I_k^{\text{ion}}, & \text{on } \Gamma \\ \vec{\nabla}_{\vec{x}} \cdot \underline{j}_r^{(\alpha)} = 0, & \text{in } \Omega^{(\alpha)}, \\ \underline{j}_r^{(\alpha)} = -\left(\sigma^{(\alpha)} \vec{\nabla}_{\vec{x}} u^{(\alpha)} + \mathrm{F} \sum_{k=1}^{N^{(\alpha)}} z_k D_k \vec{\nabla}_{\vec{x}} c_k^{(\alpha)}\right), & \text{in } \Omega^{(\alpha)}, \\ s_{\Gamma} \underline{j}_r^{(i)} \cdot \underline{n}^{(i)} = -s_{\Gamma} \underline{j}_r^{(e)} \cdot \underline{n}^{(e)} = C_m \partial_t u^{(m)}|_{\xi} + I^{\text{ion}}, & \text{on } \Gamma, \end{array} \right. \quad (1.5) $$

where $s_{\Gamma} = J \|F^{-T} \cdot \underline{n}_0^{(\alpha)}\|$ represents the stretch ratio of the membrane in the current configuration with respect to the reference one, $u^{(m)} = u^{(i)} - u^{(e)}$ denotes the transmembrane potential and $C_m$ the capacitance per unit area of the membrane, see [11] for all details. In addition, the auxiliary unknowns $\mu_k^{(\alpha)}$ and $\mu^{(\alpha)}$ are given by $\mu_k^{(\alpha)} = (z_k c_k^{(\alpha)} \mu^{(\alpha)})/(\sum_{\ell} (z_{\ell})^2 c_{\ell}^{(\alpha)} \mathrm{F})$ and $\mu^{(i)} = -\mu^{(e)} = C_m u^{(m)}$ on $\Gamma$.

### 1.1.4 Bidomain model simplifications

In order to derive the classical bidomain equations, it is assumed that diffusion can be neglected compared to drift in ionic fluxes [Mor09]

$$ \left\{ \begin{array}{ll} \vec{\nabla}_{\vec{x}} \cdot (\sigma^{(\alpha)} \vec{\nabla}_{\vec{x}} u^{(\alpha)}) = 0, & \text{in } \Omega^{(\alpha)} \\ \sigma^{(i)} \vec{\nabla}_{\vec{x}} u^{(i)} \cdot \underline{n}^{(i)} = \sigma^{(e)} \vec{\nabla}_{\vec{x}} u^{(e)} \cdot \underline{n}^{(i)}, & \text{on } \Gamma \\ -s_{\Gamma} \sigma^{(i)} \vec{\nabla}_{\vec{x}} u^{(i)} \cdot \underline{n}^{(i)} = C_m \partial_t u^{(m)}|_{\xi} + I^{\text{ion}}, & \text{on } \Gamma. \end{array} \right. \quad (1.6) $$

Whether or not this assumption is valid in practice is an open question that goes beyond the scope of this work. Finally, one can observe that the previous equations define $u^{(i)}$ or $u^{(e)}$ up to the same constant. Therefore, we choose to impose $\int_{\Gamma} u^{(e)} d\gamma = 0$. To conclude the description of the microscopic bidomain model, the last subsection is dedicated to the description of the transmembrane currents.
### 1.1.5 Ionic models and mechanical effects

**Physiological or phenomenological ionic models** Physiological ionic models aim at describing ionic exchanges up to some varying degree of detail. The first such model was proposed for neurons by Hodgkin and Huxley [HH52]. Some models were later introduced to more specifically address cardiac electrophysiology with – to cite just a few classical ones – the models proposed in [CRN98, NFF+98] for atrial cells and in [tTNNP04, tTP06, OVVR11] for ventricular cells. By contrast, so-called phenomenological ionic models are intended to focus on reproducing the variations of the transmembrane potential with much fewer state variables. The first one was independently proposed by FitzHugh and Nagumo [Fit61, NAY62] for neurons, and was later adapted to cardiac cells in [BCP94, AP96]. More recent phenomenological model were proposed in [MS03, BOCF08]. Since most ionic models can be rewritten as follows

$$I^{\text{ion}} := I^{\text{ion}}(u^{(m)}, w), \quad \partial_t w + g(u^{(m)}, w) = 0,$$

where $w$ corresponds to one (or more) additional variable(s) called state variable(s), we will consider this formulation in what follows, see [8] for more details.

**Stretch-Activated Channels** To a first approximation, cardiac electrophysiology and mechanics are often assumed to be weakly coupled, that is, electrophysiology is considered to be an input that drives the active part of mechanical behavior. However, it is known that mechanical deformation alters electrical properties. This is the so-called mechano-electrical feedback (MEF), which is thought to be involved in cardiac arrhythmias and other pathologies, see, e.g., [ABE16] and references therein. One of the mechanisms of MEF is provided by the stretch-activated channels (SAC). To model this effect, some parts of ionic models are made to depend on deformation, which usually leads to a decrease in action potential duration and amplitude and an increase in resting transmembrane potential [HM05, AHZ+13, PKN05, AHZ+13, AHZ+13]. However, as the role and mechanisms of SAC are still not fully understood [WGRS12], we disregard the direct dependence on mechanical deformations of ionic model to evaluate other effects.

## 1.2 Mathematical analysis of the periodic homogenization

In this section we neglect the mechanical deformations, which means that we consider System (1.6) in which $\phi$ is the identity, which means that $\underline{x}(t) = \xi, \forall t, s_{\Gamma} = 1$, and $\underline{v}_s = 0$.

This section is a summary of our article [8]. The goal of our work is to fill a gap in the literature by providing a complete mathematical analysis based on the 2-scale convergence theory of the homogenization process leading to the macroscopic bidomain model. Our analysis is exhaustive in the sense that we provide existence and uniqueness results, non-dimensionalization of equations, and 2-scale convergence results – in particular, for the nonlinear terms supported on the membrane surface – within the same mathematical framework. To anticipate reasonable modeling assumptions, we assume that the electric conductivities are tensorial and spatially variable at the microscopic scale. We also consider ionic models of both types (physiological and phenomenological) that can vary uniformly in space (to consider ventricular or atrial cells, for example). We carefully introduce the various standard assumptions satisfied by ionic terms and distinguish the models that are compatible with our analysis. We believe that this work will enable the analysis of more complex models by laying the foundation for the 2-scale analysis of the
bidomain equations. More specifically, we study the following system

$$ 
\left\{
\begin{array}{l}
\displaystyle \vec{\nabla}_{\vec{x}} \cdot (\vec{\sigma}^{(\alpha)} \vec{\nabla}_{\vec{x}} u^{(\alpha)}) = 0 \qquad \Omega^{(\alpha)}, \\[1.5em]
\displaystyle \vec{\sigma}^{(i)} \vec{\nabla}_{\vec{x}} u^{(i)} \cdot \vec{n}^{(i)} = \vec{\sigma}^{(e)} \vec{\nabla}_{\vec{x}} u^{(e)} \cdot \vec{n}^{(i)} \qquad \Gamma^{(m)}, \\[1.5em]
\displaystyle -\vec{\sigma}^{(i)} \vec{\nabla}_{\vec{x}} u^{(i)} \cdot \vec{n}^{(i)} = C_m u^{(m)} + I^{\text{ion}}(u^{(m)}, w) \qquad \Gamma^{(m)}, \\[1.5em]
\displaystyle u^{(m)} = u^{(i)} - u^{(e)} \qquad \Gamma^{(m)}, \\[1.5em]
\displaystyle \partial_t w = -g(u^{(m)}, w) \qquad \Gamma^{(m)}, \\[1.5em]
\displaystyle \vec{\sigma}^{(\alpha)} \vec{\nabla}_{\vec{x}} u^{(e)} \cdot \vec{n}^{(\alpha)} = 0 \qquad \partial\Omega^{(\alpha)} \cap \partial\Omega, \\[1.5em]
\displaystyle \int_{\Gamma^{(m)}} u^{(e)} d\gamma = 0.
\end{array}
\right.
\tag{1.7}
$$

Among the modeling ingredients that might fit our context, one might consider: heterogeneous concentrations of ionic species within cells, gap junctions [HP10], and microscopic cardiac fiber structure in the context of local 2-scale convergence [Bri93, Pta17]. In Section 1.2.1, the main results of the existence and uniqueness of the heterogeneous microscopic bidomain model are given. In Section 1.2.2 the main results of the homogenization process of the heterogeneous microscopic bidomain model are given.

### 1.2.1 Analysis of the microscopic bidomain model

In [8], the assumptions used to prove existence and uniqueness of the microscopic bidomain model for various ionic models are very carefully selected.

**Assumption 1.** Summary of the assumptions:

(a) *Regularity of the diffusion tensors.* $\vec{\sigma}^{(\alpha)}$ are bounded, symmetric, definite, positive and coercive and allow to define norms in $L^2(\Omega^{(\alpha)})$.

(b) *Regularity of the initial conditions.* $\overline{u^{(m)}}^0 \in L^2(\Gamma^{(m)})$ and $\overline{w}^0 \in L^2(\Gamma^{(m)})$.

(c) *Regularity condition of the ionic terms.* $I^{\text{ion}} \in C^0(\bar{\Omega} \times \mathbb{R}^2)$ and $g \in C^0(\bar{\Omega} \times \mathbb{R}^2)$.

(d) *Growth conditions of the ionic terms.* $\exists C_\infty > 0, \forall \vec{x} \in \Omega, \forall (v, w) \in \mathbb{R}^2$,

$$ |I^{\text{ion}}(\vec{x}, v, w)| \leq C_{\infty}(|v|^3 + |w| + 1) \quad \text{and} \quad |g(\vec{x}, v, w)| \leq C_{\infty}(|v|^2 + |w| + 1). \tag{1.8} $$

(e) *Signed condition of the ionic terms.* There exist $\mu > 0$ and $C_I > 0$ such that for all $\vec{x} \in \Omega$ and $(v, w) \in \mathbb{R}^2$, we have

$$ v I^{\text{ion}}(\vec{x}, v, w) + \mu w g(\vec{x}, v, w) \geq -C_I (|v|^2 + |w|^2 + 1). \tag{1.9} $$

(f) *Compensation of the lack of regularity on the gating variables.* Depending on whether the model is physiological or phenomenological, one of the following assumptions applies

(i) *Global lipschitz property.* There exists a positive scalar $L_g > 0$, such that for all $\vec{x} \in \Omega$, $(v_1, v_2) \in \mathbb{R}^2$ and $(w_1, w_2) \in \mathbb{R}^2$,

$$ |g(\vec{x}, v_1, w_1) - g(\vec{x}, v_2, w_2)| \leq L_g |v_1 - v_2| + L_g |w_1 - w_2|. \tag{1.10} $$
(ii) Decomposition of the nonlinear terms. There exist continuous functions

$$ (f_1, f_2, g_1, g_2) \in [C^0(\Omega \times \mathbb{R}^2)]^4 $$

such that

$$ \begin{cases} I^{\text{ion}}(\vec{x}, v, w) = f_1(\vec{x}, v) + f_2(\vec{x}, v) w, \\ g(\vec{x}, v, w) = g_1(\vec{x}, v) + g_2(\vec{x}) w, \end{cases} $$

and there exist positive constants $C_1$, $c_1$ and $C_2$ such that for all $\vec{x} \in \Omega$ and $v \in \mathbb{R}$,

$$ v f_1(\vec{x}, v) \geq C_1 |v|^4 - c_1 (|v|^2 + 1) \quad \text{and} \quad |f_2(\vec{x}, v)| \leq C_2 (|v| + 1). \quad (1.11) $$

(g) One-sided Lipschitz condition for the uniqueness of a solution. There exist $\mu > 0$ and $L_I > 0$ such that for all $\vec{x} \in \Omega$, $(v_1, w_1) \in \mathbb{R}^2$ and $(v_2, w_2) \in \mathbb{R}^2$

$$ \begin{aligned} & (I^{\text{ion}}(\vec{x}, v_1, w_1) - I^{\text{ion}}(\vec{x}, v_2, w_2))(v_1 - v_2) + \mu (g(\vec{x}, v_1, w_1) - g(\vec{x}, v_2, w_2))(w_1 - w_2) \\ & \geq -L_I (|v_1 - v_2|^2 + |w_1 - w_2|^2). \end{aligned} \quad (1.12) $$

Assumptions from 1(a) to 1(f) are satisfied for classical ionic models, see Table 1.1 for a classification. Assumption 1(g) is not satisfied for the Aliev-Panfilov [AP96], the MacCulloch [BCP94], and the Mitchell-Schaeffer [MS03] models. The uniqueness of a solution for these models is still an open problem.

<table><caption>Table 1.1 – Ionic models and verified assumptions</caption><thead><tr><th>Assumption</th><th>1(d)</th><th>1(e)</th><th>1(f)</th><th>1(g)</th></tr></thead><tbody><tr><td>FitzHugh-Nagumo</td><td>✓</td><td>✓</td><td>(i)</td><td>✓</td></tr><tr><td>Roger MacCulloch</td><td>✓</td><td>✓</td><td>(i)</td><td></td></tr><tr><td>Aliev-Panfilov</td><td>✓</td><td>✓</td><td>(ii)</td><td></td></tr><tr><td>Regularized Mitchell-Schaeffer</td><td>✓</td><td>✓</td><td>(i)</td><td></td></tr><tr><td>Physiological models</td><td>✓</td><td>✓</td><td>(i)</td><td>✓</td></tr></tbody></table>

### 1.2.1.1 Existence and uniqueness analysis

In the literature, the analysis of the classical bidomain model is mostly performed at the macroscopic scale, see [CFPS14] for a presentation of the different approaches. In any case, to justify the homogenization process, a mathematical analysis at the microscopic scale is required. In this section, we give existence and uniqueness results for solutions of System (1.7) using the Faedo-Galerkin approach. The main five steps are quickly detailed below.

**Step 1 - Elimination of the unknown electric potentials** Following an idea developed in [BCP09] for the bidomain equation at the macroscopic level or in [KLPW14] for an electroporation model at the microscopic level, one can prove, see [8, Section 3.1], that there exists a regular linear operator $\mathcal{A}$ defined as a combination of Dirichlet-to-Neumann type of operators such that solutions of System (1.7) verify

$$ \left\{ \begin{array}{l} C_m \partial_t u^{(m)} + I^{\text{ion}}(u^{(m)}, w) = -\mathcal{A}(u^{(m)}), \quad \Gamma^{(m)}, \\ \partial_t w = -g(u^{(m)}, w), \quad \Gamma^{(m)}. \end{array} \right. \quad (1.13) $$
**Step 2 - Evolution equation of the gating variable.** Let $U \in C^0([0, T]; L^2(\Gamma^{(m)}))$, a given electrical potential. The term $g(u^{(m)}, w)$ in (1.13) is replaced by $g(U, w)$ with $w = w_U$ the corresponding solution.

**Lemma 1.** If Assumption 1(b), 1(c), 1(d) and 1(f)(i) hold, there exists a unique function $w_U \in C^1([0, T]; L^2(\Gamma^{(m)}))$, which is a solution of

$$ \left\{ \begin{array}{ll} \partial_t w_U + g(\vec{x}, U, w_U) = 0, & \Gamma^{(m)}, \forall t \in [0, T], \\ w_U(\vec{x}, 0) = \bar{w}^0(\vec{x}), & \Gamma^{(m)}. \end{array} \right. \qquad (1.14) $$

and there exists a positive scalar $C_w$, that depends only on $T$ and $L_g$, such that for any function $V \in C^0((0, T); L^2(\Gamma^{(m)}))$,

$$ \int_{0}^{T} \|w_{U}(t) - w_{V}(t)\|_{L^{2}(\Gamma^{(m)})}^{2} \leq C_{w} \int_{0}^{T} \|U(t) - V(t)\|_{L^{2}(\Gamma^{(m)})}^{2} dt. $$

The proof of Lemma 1 is done by considering smooth approximations of $U$ and $\bar{w}^0$. Then, the problem reduces to the analysis of an ODE where the space variable plays the role of a parameter. Finally, the solution of (1.14) is constructed by a limit process using the density of smooth functions into $L^2(\Gamma^{(m)})$.

### Step 3 - Existence result for the microscopic bidomain equation

**Theorem 1.** If Assumptions from 1(a) to 1(f) hold, there exist

$$ u^{(m)} \in C^0([0, T]; L^2(\Gamma^{(m)})) \cap L^2((0, T); H^{1/2}(\Gamma^{(m)})), \quad \partial_t u^{(m)} \in L^2((0, T); H^{-1/2}(\Gamma^{(m)})), $$

and $w \in H^1((0, T); L^2(\Gamma^{(m)}))$, which are solutions of

$$ \left\{ \begin{array}{ll} C_m \partial_t u^{(m)} + A(u^{(m)}) + I^{ion}(u^{(m)}, w) = 0, & \tilde{H}^{-1/2}(\Gamma^{(m)}), \quad \text{a.e. } t \in (0, T), \\ \partial_t w + g(u^{(m)}, w) = 0, & \Gamma^{(m)}, \quad \text{a.e. } t \in (0, T), \end{array} \right. \qquad (1.15) $$

$$ \text{and } u^{(m)}(\vec{x}, 0) = \overline{u^{(m)}}^{0}(\vec{x}), w(\vec{x}, 0) = \bar{w}^{0}(\vec{x}), \Gamma^{(m)}. $$

The proof of Theorem 1 is done using the Faedo-Galerkin method. More precisely, the equations are first space-discretized using a finite dimensional basis of $L^2(\Gamma^{(m)})$ constructed with the eigenvectors of $A$. After the discretization, it is proven that semi-discrete solutions exist by applying the Cauchy-Peano theorem (to be more specific, we use the more general Carathéodory's existence theorem) on systems of ordinary differential equations. Finally, by a limit procedure, the existence of solutions is proven for the weak form of (1.15).

### Step 4 - Uniqueness results for the microscopic bidomain equation
Uniqueness is proven by standard energy techniques for models satisfying Assumption 1(g).

### Step 5 - Post-processing of the intra- and extra-cellular potentials
From the solution $u^{(m)} = u^{(i)} - u^{(e)}$ given by Theorem 1, we can recover the intra- and extra-cellular potentials

$$ u^{(i)} \in L^2((0, T); H^1(\Omega_0^{(i)})) \text{ and } u^{(e)} \in L^2((0, T); H^1(\Omega_0^{(e)})) $$

following the strategy given in Section 3.4 of [8].
### 1.2.2 Homogenization of the bidomain equations

**Context** After existence and uniqueness, we are now in a position to perform a rigorous homogenization of the microscopic bidomain model (1.6). The microscopic model is useless for the whole heart in terms of numerical applications. At the macroscopic scale, the heart appears as a continuous material with a fiber-based structure. At this scale, the intracellular and extracellular media are indistinguishable. Our goal is to use a homogenization of the microscopic bidomain model to obtain a bidomain model in which all unknowns are defined everywhere. Formally, the heart volume after homogenization is “$\Omega = \Omega^{(i)} = \Omega^{(e)}$” (for simplicity, the subscript 0 has been removed but as the mechanical deformations are not considered, it corresponds to the reference domain). Homogenization of partial differential systems is a well-known technique, see [BLP78] for a reference book on this subject. Here the medium is considered periodic of period $\epsilon$, and equations are derived by asymptotic analysis with respect to $\epsilon$. The periodicity allows to use the results of 2-scale convergence, see [All92], a method which has been used in many areas of science and engineering. We denote by $Y$ the open reference domain that is used to define the idealized micro-structure corresponding to the periodic arrangement of cardiac cells. This micro-structure is decomposed into two open connected subdomains: the intracellular part $Y^{(i)}$ and the extracellular part $Y^{(e)}$. We have $Y_i \cap Y_e = \emptyset$, $\bar{Y} = \bar{Y}_i \cup \bar{Y}_e$. The intra- and the extra-cellular domains are separated by $\Gamma_Y$. The global position vector is denoted by $\vec{x}$ and the local position vector by $\vec{t}$. In this context, the domain $\bar{\Omega} = \bar{\Omega}_{\epsilon}^{(i)} \cup \bar{\Omega}_{\epsilon}^{(e)}$ is then defined by $\epsilon$-periodicity and we denote by $\Gamma_{\epsilon}^{(m)}$ the boundary between the intra- and the extra-cellular domains $\Omega_{\epsilon}^{(i)}$ and $\Omega_{\epsilon}^{(e)}$ i.e. $\Omega_{\epsilon}^{(\alpha)} = \cup_k (\epsilon Y^{(\alpha)} + \epsilon \vec{w}_k)$ and then $\Gamma_{\epsilon}^{(m)} = \cup_k (\epsilon \Gamma_Y + \epsilon \vec{w}_k)$, where $\vec{w}_k$ (with $\vec{w}_0 = \vec{0}$) is the vector corresponding to the translation between the considered cell and the reference cell. Note that by construction, we have $\vec{t} = \vec{x}/\epsilon$. We assume that the diffusion tensors depend on the two scales i.e. $\vec{\sigma}_{\epsilon}^{(\alpha)}(\vec{x}) = \vec{\sigma}^{(\alpha)}(\vec{x}, \frac{\vec{x}}{\epsilon})$.

A classical article for the formal homogenization of the microscopic bidomain model – when the conductivities $\sigma^{(\alpha)}$ are strictly positive constants – is [NK93]. In [CFS02], the homogenization of the microscopic bidomain equations (with constants conductivities) is presented using formal asymptotic analysis. It should also be noted that in [ACFS00, PSCF05], the same type of results have been proven using the theory of $\Gamma$-convergence in some simplified situations. The approach presented below uses the 2-scale convergence method, see [All92], to extend the results obtained in [ACFS00, PSCF05]. As typical for this kind of homogenization problem, we adopt the following approach, see [8]:

**Step 1 - Nondimensionalization of the problem.** The microscopic bidomain equations are scaled in space and time, and written unitless. Using the characteristic values of the physical parameters of our problem (cells size, conductivities, ionic current,...), a small parameter $\epsilon$ is introduced in the equations and in the geometry. We obtain in [8, Section 4.1] the following system

$$ 
\left\{
\begin{array}{l@{\qquad}l}
\vec{\nabla}_{\vec{x}} \cdot (\vec{\sigma}_{\epsilon}^{(\alpha)} \vec{\nabla}_{\vec{x}} u_{\epsilon}^{(\alpha)}) & = 0 & \Omega_{\epsilon}^{(\alpha)}, \\
\vec{\sigma}_{\epsilon}^{(i)} \vec{\nabla}_{\vec{x}} u_{\epsilon}^{(i)} \cdot \vec{n}^{(i)} & = \vec{\sigma}_{\epsilon}^{(e)} \vec{\nabla}_{\vec{x}} u_{\epsilon}^{(e)} \cdot \vec{n}^{(i)}, & \Gamma_{\epsilon}^{(m)} \\
\vec{\sigma}_{\epsilon}^{(i)} \vec{\nabla}_{\vec{x}} u_{\epsilon}^{(i)} \cdot \vec{n}^{(i)} & = -\epsilon C_m \partial_t (u_{\epsilon}^{(i)} - u_{\epsilon}^{(e)}) - \epsilon I^{\text{ion}}(u_{\epsilon}^{(i)} - u_{\epsilon}^{(e)}, w_{\epsilon}) & \Gamma_{\epsilon}^{(m)}, \\
\partial_t w_{\epsilon} & = -g(u_{\epsilon}^{(i)} - u_{\epsilon}^{(e)}, w_{\epsilon}), & \Gamma_{\epsilon}^{(m)}. \\
\vec{\sigma}_{\epsilon}^{(i)} \vec{\nabla}_{\vec{x}} u_{\epsilon}^{(\alpha)} \cdot \vec{n}^{(\alpha)} & = 0, & \partial\Omega \cap \partial\Omega_{\epsilon}^{(\alpha)}, \\
u_{\epsilon}^{(i)}(\cdot, 0) - u_{\epsilon}^{(e)}(\cdot, 0) & = \overline{V_{m,\epsilon}}^0, & \Gamma_{\epsilon}^{(m)} \\
w_{\epsilon}(\cdot, 0) & = \overline{w}_{\epsilon}^0, & \Gamma_{\epsilon}^{(m)}, \\
\int_{\Gamma_{\epsilon}^{(m)}} u_{\epsilon}^{(e)} d\gamma & = 0.
\end{array}
\right.
$$
**Step 2 - Uniform estimate of solutions.** Using energy estimates, we prove that the norms of the solutions – as well as norms of the nonlinear terms – are uniformly bounded with respect to $\varepsilon$ if the initial conditions follow uniform estimates, see [8, Section 4.2].

**Step 3 - Two-scale convergence.** The limit equations are deduced by application of the 2-scale convergence theory. One of the main difficulties of this step is the convergence analysis of the nonlinear terms. We quickly detail this part below. In [8], one can find a detailed proof of the 2-scale convergence of $u_{\varepsilon}^{(\alpha)}$ to the sought macro-scale solutions $u_{0}^{(\alpha)}$ verifying

$$ 
\left\{ 
\begin{array}{l}
-\nabla_{\vec{x}} \cdot \vec{T}^{(i)} \nabla_{\vec{x}} u_{0}^{(i)} + A_{m} C_{m} (u_{0}^{(i)} - u_{0}^{(e)}) + A_{m} I^{\text{ion}}(u_{0}^{(i)} - u_{0}^{(e)}, w_{0}) = 0, \quad \Omega \times (0, T), \\
-\nabla_{\vec{x}} \cdot \vec{T}^{(e)} \nabla_{\vec{x}} u_{0}^{(e)} - A_{m} C_{m} (u_{0}^{(i)} - u_{0}^{(e)}) - A_{m} I^{\text{ion}}(u_{0}^{(i)} - u_{0}^{(e)}, w_{0}) = 0, \quad \Omega \times (0, T), \\
\partial_{t} w_{0} + g(u_{0}^{(i)} - u_{0}^{(e)}, w_{0}) = 0, \quad \Omega \times (0, T), \\
(\vec{T}^{(\alpha)} \cdot \nabla_{\vec{x}} u_{0}^{(\alpha)}) \cdot \vec{n} = 0, \quad \partial\Omega \times (0, T), \\
(u_{0}^{(i)} - u_{0}^{(e)})(\vec{x}, 0) = \overline{u^{(m)}(\vec{x})}, \quad \Omega, \\
w_{0}(\vec{x}, 0) = \overline{w^{0}(\vec{x})}, \quad \Omega, \\
\int_{\Omega} u_{0}^{(e)} d\vec{x} = 0
\end{array}
\right. 
$$

The tensors $\vec{T}^{(\alpha)}$, $\alpha = (i, e)$ are the associated effective medium tensors

$$ (\vec{T}^{(\alpha)})_{j,k} = \frac{1}{|Y|} \int_{Y^{(\alpha)}} \vec{\sigma}^{(\alpha)} (\vec{\nabla}_{\vec{y}} \mathcal{X}_{j}^{(\alpha)} + \vec{e}_{j}) \cdot (\vec{\nabla}_{\vec{y}} \mathcal{X}_{k}^{(\alpha)} + \vec{e}_{k}) d\vec{t}, \quad (1.16) $$

where the canonical functions $\mathcal{X}_{j}^{(\alpha)}$, $j = 1..3$ belong to $H_{\#}^{1}(Y^{(\alpha)})$ and are uniquely defined by the following variational formulation – called the cell problem

$$ 
\left\{ 
\begin{array}{l}
(\vec{\sigma}^{(\alpha)}(\vec{e}_j + \vec{\nabla}_{\vec{y}} \mathcal{X}_j^{(\alpha)}), \vec{\nabla}_{\vec{y}} \psi)_{Y^{(\alpha)}} = 0, \quad \forall \psi \in H_{\#}^{1}(Y^{(\alpha)}), \\
\int_{Y^{(\alpha)}} \mathcal{X}_j^{(\alpha)} d\vec{t} = 0.
\end{array}
\right. 
\quad (1.17) $$

For the proof, we use the 2-scale convergence theory developed in [All92] and [Ngu89]. The treatment of the nonlinear ionic terms represents one of the main difficulties. To solve this problem, we present an approach using the general ideas of [All92]. It is based on the one-sided Lipschitz assumption (Assumption 1(g)). The final difficulty is that the nonlinear terms are at the boundary of the domain. The 2-scale convergence theory must be adapted, and for this we use the results given in [ADH95].

Since the analysis has already been performed for the microscopic bidomain model, we can infer – up to a subsequence – the existence of a solution for the macroscopic bidomain model. The proposed model is a generalization of the very classical macroscopic bidomain model in which constant electric conductivities are considered. Compared to previous studies, see for example [ACFS00] and [PSCF05], and in order to anticipate meaningful modeling assumptions, we have considered at the microscopic level that the electric conductivities are tensorial. This appears in the definition of the cell problem (1.17) and in the definition of the tensor (1.16). Moreover, we have shown that the classical macroscopic bidomain model formally obtained in [CFS02] and [NK93] is valid for the ionic terms under some more general conditions than those assumed in [ACFS00] and [PSCF05]. More precisely, we extended the validity of the macroscopic bidomain equations to spatially variable physiological models.
## 1.3 Effects of mechanical deformations on the bidomain model

### 1.3.1 Macroscopic bidomain incorporating the mechanical deformations

Now we return to the incorporation of the mechanical deformations. We use our expertise acquired in the full analysis of the modeling when the mechanical deformations are neglected to perform a formal asymptotic analysis on System (1.5), which we rewrite here in the Lagrangian framework, i.e. in terms of the spatial variable associated with the reference configuration. For this purpose, it is classical, as is usual in the derivation of conservation laws in moving domains [TM05], to introduce the following modified concentration and (relative) flux quantities

$$C_k^{(\alpha)} = J c_k^{(\alpha)}, \quad J_k^{(\alpha)} = J \underline{F}^{-1} \cdot \underline{J}_{k,r}^{(\alpha)}.$$

Our homogenization Ansatz is that each main unknown of the microscopic problem defined on the reference domain is the sum of two terms defined everywhere in $\Omega_0$:

* a leading term independent of $\varepsilon$,

* correction terms dependent on $\varepsilon$ and with local variations inside each cell.

Regarding the membrane equations, we consider asymptotic assumptions consistent with the nondimensionalized system presented in the previous section meaning that they are scaled by like $\varepsilon$.

We then incorporate the equations of our homogenization Ansatz into the governing equations and we identify the equations satisfied by the leading order terms $C_{k,0}^{(\alpha)}$ and $u_0^{(\alpha)}$. It should be noted, however, that these terms will be mixed with higher order terms in the equations that we will directly infer, due to space differentiation with the coexistence of the $\xi$ and $\xi/\varepsilon$ variables in various terms.

We only give the results but all the computations are in [11]. First, we introduce the basic cell problem in which, for any vector $\underline{v} \in \mathbb{R}^3$, we call $\mathcal{X}^{(\alpha)}(\underline{v})$ the periodic solution (defined up to a constant) of

$$\left\{ \begin{array}{ll} \vec{\nabla}_{\vec{\eta}} \cdot \left[ \left( \underline{F}^T \cdot \underline{F} \right)^{-1} \cdot \vec{\nabla}_{\vec{\eta}} \mathcal{X}^{(\alpha)}(\underline{v}) \right] = 0, & \text{in } Y^{(\alpha)} \\ \left[ \left( \underline{F}^T \cdot \underline{F} \right)^{-1} \cdot \vec{\nabla}_{\vec{\eta}} \mathcal{X}^{(\alpha)}(\underline{v}) \right] \cdot \underline{n}_Y^{(\alpha)} = - \left[ \left( \underline{F}^T \cdot \underline{F} \right)^{-1} \underline{v} \right] \cdot \underline{n}_Y^{(\alpha)}, & \text{on } \Gamma_Y \end{array} \right. \quad (1.18)$$

Using this solution and as it is classical in periodic homogenization analysis (see. [All12, Chapter 1]), we introduce the tensor $\underline{T}^{(\alpha)}$ defined by

$$\underline{T}^{(\alpha)} \cdot \underline{u} = \frac{1}{|Y^{(\alpha)}|} \int_{Y^{(\alpha)}} (\underline{F}^T \cdot \underline{F})^{-1} \cdot (\underline{u} + \vec{\nabla}_{\vec{\eta}} \mathcal{X}^{(\alpha)}(\underline{u})) d\underline{\eta}. \quad (1.19)$$

Note that $\underline{T}^{(\alpha)}$ only depends on the geometry of the cell and on the deformation through $\underline{F}$. We obtain

$$\left\{ \begin{array}{l} \phi^{(\alpha)} \partial_t C_{k,0}^{(\alpha)} + \vec{\nabla}_{\vec{\xi}} \cdot \phi^{(\alpha)} \underline{J}_{k,0}^{(\alpha)} + A_m p_{k,0}^{(\alpha)} = 0, \\ \underline{J}_{k,0}^{(\alpha)} = -J D_k \underline{T}^{(\alpha)} \cdot \left( \vec{\nabla}_{\vec{\xi}} (J^{-1} C_{k,0}^{(\alpha)}) + \frac{q z_k (J^{-1} C_{k,0}^{(\alpha)})}{k_B T} \vec{\nabla}_{\vec{\xi}} u_0^{(\alpha)} \right), \end{array} \right. \quad (1.20)$$

with $\phi^{(\alpha)} = \frac{|Y^{(\alpha)}|}{|Y|}$, the volume fractions of intra- and extra-cellular media in the reference configuration and

$$p_{k,0}^{(i)} = \mu_{k,0}^{(i)} + I_{k,0}^{\text{ion}}, \quad p_{k,0}^{(e)} = \mu_{k,0}^{(e)} - I_{k,0}^{\text{ion}}.$$
<figure><img src="image_13.png" /><figcaption>Figure 1.2 – Left: Mesh of the ventricles. Right: Fibers in the ventricles.</figcaption></figure>

The additional relations pertaining to surface charges read

$$ \mu_{k,0}^{(\alpha)} = \frac{z_k C_{k,0}^{(\alpha)}}{\sum_{\ell} (z_{\ell})^2 C_{\ell,0}^{(\alpha)}} \frac{\mu_0^{(\alpha)}}{F}, \quad \mu_0^{(i)} = -\mu_0^{(e)} = \hat{C}_m V_{m,0}. $$

Note that the terms $A_m p_{k,0}^{(\alpha)}$ in (1.20) represent the quantity of ions going from the volume to the membrane, for each species and in each domain. Finally, we complete this set of equations with

$$ \left\{ \begin{aligned} \vec{\nabla}_{\vec{\xi}} \cdot \phi^{(\alpha)} \hat{J}_{0}^{(\alpha)} + A_m p_{0}^{(\alpha)} &= 0, \\ \hat{J}_{0}^{(\alpha)} &= -J \underline{T}^{(\alpha)} \cdot \left( \sigma_{0}^{(\alpha)} \vec{\nabla}_{\vec{\xi}} u_{0}^{(\alpha)} + F \sum_{k=1}^{N^{(\alpha)}} z_{k} D_{k} \vec{\nabla}_{\vec{\xi}} (J^{-1} C_{k,0}^{(\alpha)}) \right), \end{aligned} \right. \quad (1.21) $$

where $p_0^{(i)} = -p_0^{(e)} = \hat{C}_m V_{m,0} + I_0^{\text{ion}}$, $\sigma_0^{(\alpha)} = Fq \sum_{k=1}^{N^{(\alpha)}} m_k(z_k)^2 J^{-1} C_{k,0}^{(\alpha)}$ and $I_0^{\text{ion}} = F \sum_{k=1}^{N} z_k I_{k,0}^{\text{ion}}$.

To obtain a macroscopic bidomain equation including the effects of mechanical deformations, we need the last hypothesis. Indeed, assuming that the diffusion is negligible compared to the drift and denoting $\Sigma_0^{(\alpha)} = J\sigma_0^{(\alpha)}\phi^{(\alpha)}\underline{T}^{(\alpha)}$, one can simplify System (1.21) into

$$ \left\{ \begin{aligned} \vec{\nabla}_{\vec{\xi}} \cdot (\Sigma_0^{(i)} \cdot \vec{\nabla}_{\vec{\xi}} u_0^{(i)} + \Sigma_0^{(e)} \cdot \vec{\nabla}_{\vec{\xi}} u_0^{(e)}) &= 0, \\ A_m \left[ \hat{C}_m \frac{\partial V_{m,0}}{\partial t} + I_0^{\text{ion}} \right] - \vec{\nabla}_{\vec{\xi}} \cdot (\Sigma_0^{(i)} \cdot \vec{\nabla}_{\vec{\xi}} u_0^{(i)}) &= 0, \end{aligned} \right. \quad (1.22) $$

Of course, these equations can also be rewritten in the Eulerian framework, see [11].
<figure><img src="image_14.png" /><figcaption>Figure 1.3 – Comparison between the results of the bidomain model with (left) and without (right) deformations.</figcaption></figure>

### 1.3.2 Numerical illustrations

In this section, we would like to make a comparison between a simulation in which the mechanical deformations are taken into account and a simulation in which they are neglected. The goal is to see if neglecting the mechanical effects is reasonable in practice, and which geometry is most appropriate when neglecting them.

All the main components required to run a simulation – as variational formulation, parameters values, geometry, mesh and fiber directions (see Figure 1.2), ionic model, external stimulus, space and time discretizations and mechanical displacements – are given in [11].

**Simulation results and discussion** We display the simulation results in Figure 1.3 and in Figure 1.4 for a cut of the heart, where we compare for each time considered the potential maps obtained with and without taking the deformations into account, albeit applying the same deformation in the visualization to ease the comparison. The simulation starts when the atria contract and eject the blood in the ventricles between $t = 0$ and 190 ms, i.e. when the volumes of the ventricles are near their maximum. The applied current $I^{app}$ stimulates the endocardium at $t = 190$ ms. As we can see in Fig. 1.3 and 1.4, the depolarization is faster when deformations are considered, which is explained by the fact that the signal mainly propagates through the ventricular wall thickness from the endocardium to the epicardium. Indeed, the wall thickness is thinner when the walls are stretched. The depolarizations are complete at 211 ms vs. 214 ms with and without deformations, respectively, which gives a delay of 3 ms on a depolarization duration of 21 ms, i.e. a difference of 15%. Near the end of depolarization (around $t = 210$ ms), ejection begins. Before the end of ejection (around $t = 365$ ms), repolarization occurs. We can see ($t = 390$ ms) that the repolarization appears somewhat faster when the deformations are considered, but the advance of 3 ms due to depolarization has to be considered. The repolarizations are complete at 410 ms vs. 411 ms with and without deformations, respectively, which gives a small difference in the end.

We can summarize our observations by saying that the effects of deformations as represented by our model are predominantly seen in the depolarization phase, where they speed up the propagation through the wall, with an impact on the depolarization timings that is comparable – in relative value – to
<figure><img src="image_15.png" /><figcaption>Figure 1.4 – Comparison between the results of the bidomain model with (left) and without (right) deformations in a cut of the heart.</figcaption></figure>

the deformation amplitude. We point out that this is also comparable to the impact obtained when only considering changes in the constitutive laws to account for the MEF. Note that changing the geometry on which the fixed-domain electrophysiology problem is solved to one closer to end-diastolic configuration would most probably reduce the discrepancy in the depolarization phase with the solution of the model incorporating mechanical deformations. Nevertheless, this would not be satisfactory in our context where only the mechanical reference configuration (close to end-systole) is assumed to be given, and moreover in this case it would instead introduce discrepancies in the repolarization phase. In other words, if conductivity parameters are calibrated in one given configuration (end-systolic or end-diastolic), using the same parameters in a phase where the deformed geometry is far from this configuration will introduce discrepancies between the solutions of the electrophysiology problem with and without deformations.

**Electrocardiogram simulations** An ECG is a recording of the electrical activity of the heart [MP95, War75] measured at the outer surface of the body. This procedure – which is non-invasive and inexpensive – is the most standard diagnosis tool for cardiologists to detect cardiac pathologies. The literature on the simulation of ECG is abundant, see e.g. [BP11, Hui98, LBG$^{+03}$, KSW$^{+07}$]. The simulation of 12-lead ECGs based on partial differential equations (PDEs) – as proposed here – appeared during the last decade [BCF$^{+10}$, MDFG12, PDG03, PDV09, TDP$^{+04}$], more recently than simulations based on cellular automata [WOH$^{+95}$].

To model an ECG, the coupling between the heart and the body must be taken into account. It can be done through a diffusion problem with well-chosen transmission conditions at the heart-body interface. Our modeling choice has been detailed in [11].

Because the ECG is the most commonly used examination to assess the heart, it is important to understand the effects of mechanical deformation on the ECG in order to assimilate this type of data to estimate for example parameters of bidomain model.

Figure 1.5 displays the simulated ECGs with (blue) and without (red) deformations. The simulated ECGs can be validated against numerous criteria used and those considered in this work are given in [11] and allows to validate our modeling. The quality if these ECGs is discussed in [6] and in [11].

The main difference between the two simulations is the duration of the R-wave in the limb leads,
<figure><img src="image_16.png" /><figcaption>Figure 1.5 – Simulated ECGs with (blue) and without (red) deformations.</figcaption></figure>

which is smaller when the displacements are considered. This is due to the faster depolarization.

Our simulations reveal a significant – although not predominant – impact of mechanical deformations on the ECGs. One natural perspective would be to compare our mechanically deformed bidomain model – which takes into account the heart deformation in conservation laws – with the classical mechano-electrical feedback model in a realistic case of ECG simulation.

## 1.4 Contributions to data assimilation aspects

As just seen, we are now able to obtain very realistic simulations of the electrophysiology of the heart. We have also seen that in most cases, the mechanical deformations can be neglected. As for all natural and physical systems, a great difficulty consists now in dealing with the many uncertain quantities which must be preserved for running model simulations for a given patient. These quantities include initial conditions and physical parameters of the model – as for example the conductivity tensors – which are difficult to measure. Fortunately, other types of measurement exist and the objective is to use these available data in order to circumvent the uncertainties corresponding to the dynamical model.

Here, we consider the propagation maps corresponding to the position of the front at different times. These maps give the position of the depolarization fronts over time. Mathematically, we assume that the data correspond to the position of the front of the transmembrane potential $u^{(m)} = u^{(i)} - u^{(e)}$. The proposed methods can be used in other application fields where models (reaction-diffusion) and data (front position) are very similar, as for fire propagation or tumor growth. Our methods are based on sequential observers, see Section B.3 of the introduction chapter for a quick presentation.

We will first start by presenting the *Luenberger* observer for reaction-diffusion models what we have
proposed in my PhD in Section 1.4.1. In Section 1.4.2, we present a state observer based on the topological gradient that we have developed afterwards. This observer is very efficient to deal with complex patterns as one could observe when dealing with real data. In Section 1.4.3, the work that we have done in collaboration with the Carmen team, a Inria team localized at Bordeaux and specialized in modeling of cardiac electrophysiology. Our aim is to evaluate the sensitivity of our estimation strategy based on synthetic data generated with a validated realistic atrial model and then inverted with simpler modeling components. In particular, the effects of defining the muscle fibers and the corresponding anisotropic conductance parameters are examined. Finally, an application of the method to real data is presented, showing promising results. In particular, the resulting simulation smooths and completes activation maps that were spatially incomplete.

This section corresponds to very applied works. It is indeed very important for me that the innovative mathematical tools I develop are actually applicable in practice. Therefore, I try to go as far as possible to the proof of concept and also to communicate about our work in the application discipline. In this objective, the presented works have been published in two proceedings [5, 12] of the FIMH conference that aims to integrate the state of the art research and novel development efforts among others in heart modeling.

### 1.4.1 Context: a state observer using shape derivatives

**Presentation** During my PhD, we proposed, mathematically validated and numerically illustrated with synthetic data a Luenberger observer for reaction-diffusion models in the case where the uncertainties are reduced to the initial condition, see [4] for a full presentation. Following the definitions of Section B.3 of the introduction chapter, we consider the following target model

$$ \left\{ \begin{array}{l} \partial_t u = \nabla \cdot (\sigma \cdot \nabla u) + k f_r(u), \quad [0, T] \times \Omega, \\ u(0, x) = u^0(x) + \zeta_u(x), \quad \Omega, \end{array} \right. $$

where $f_r$ is a reaction term, $k$ a modeling coefficient and $u^0$ (resp. $\zeta_u$) the known (resp. the unknown) part of the initial condition. We assume that the problem is well posed i.e. there exists a unique solution $u \in L^2(0, T; \mathcal{X})$, where $\mathcal{X}$ is for example $H_0^1(\Omega)$ if we consider homogenous Dirichlet conditions. The data correspond to the position of the front of $u$ meaning that the observer operator is

$$ h: \begin{array}{l} L^2(0, T; \mathcal{X}) \\ u \end{array} \to \begin{array}{l} L^2(0, T; \mathcal{X}), \\ \phi_u, \end{array} $$

where $\phi_u$ is a level set associated with $u$ and a certain threshold value $c_{th}$ i.e. $\phi_u(t, x) \ge 0$ if and only if $u(t, x) \ge c_{th}$. One possible choice for $\phi_u$ is $(t, x) \mapsto u(t, x) - c_{th}$ but we will show in what follows that the front in the data needs to be sufficient sharp to have the convergence. Assuming that there is no measurement error, the observations become $y = h(u)$. In [4], we propose the following observer model

$$ \left\{ \begin{array}{l} \partial_t \hat{u} = \nabla \cdot (\sigma \cdot \nabla \hat{u}) + k f_r(\hat{u}) + \gamma g_u(y, \hat{u}), \quad [0, T] \times \Omega, \\ \hat{u}(0, x) = u^0(x), \quad \Omega. \end{array} \right. $$

where the gain observer is defined by

$$ g_u(y, \hat{u}) = \delta_{\{x \in \Omega, \hat{u}(t,x)=c_{th}\}}(t) \alpha(\nabla \hat{u}) ((y - C_{max}(y, \hat{u}))^2 - (y - C_{min}(y, \hat{u}))^2), $$
where

* $\delta_{\{x \in \Omega, \tilde{u}(t,x) = c_{th}\}}$ is a Dirac delta-function that ensures that the correction is only active along the front to correct,

* $\alpha$ is a strictly positive function that could be optimized to accelerate the convergence,

* and $C_{min} = \min(C_{in}, C_{out})$ and $C_{max} = \max(C_{in}, C_{out})$ with

$$C_{in}(y, \hat{u}) = \frac{\int_{\hat{u}>c_{th}} y(t, x) dx}{\int_{\hat{u}<c_{th}} dx} \quad \text{and} \quad C_{out}(y, \hat{u}) = \frac{\int_{\hat{u}<c_{th}} y(t, x) dx}{\int_{\hat{u}<c_{th}} dx}.$$

One could remark that $C_{in}$ measures the level of agreement between the observed and simulated interior areas and $C_{out}$ measures the level of disagreement between the observed and simulated exterior areas. To prove the convergence, we study the following error system

$$\begin{cases} \partial_t \tilde{u} = \nabla \cdot (\sigma \cdot \nabla \tilde{u}) + k(f_r(u) - f_r(\hat{u})) - \gamma g_u(y, \hat{u}), & [0, T] \times \Omega, \\ \tilde{u}(0, x) = \zeta_u(x), & \Omega, \end{cases}$$

whose the energy estimate is

$$\frac{d}{dt} \left( \frac{1}{2} \|\tilde{u}\|_{L^2(\Omega)}^2 \right) = - \int_{\Omega} \sigma \cdot \nabla \tilde{u} \cdot \nabla \tilde{u} \, dx + \int_{\Omega} (k(f_r(u) - f_r(\hat{u}))) \tilde{u} \, dx - \gamma \int_{\Omega} g_u(y, \hat{u}) \tilde{u} \, dx.$$

In a general context, one would like to prove the exponential or linear decrease of the energy. However, a complete mathematical analysis of the observer error is here out of the scope, due to the complexity induced by both the nonlinear reaction term and the observer correction. In [4], we were satisfied with proving the negativity of the linearized version of the observer term proving that the gain term has a stabilization effect.

**Theorem 2.** Under a condition of a sufficient sharpness of the data front detailed in [4] and when the measurement error $e$ is identically zero, the observer term is energy-decreasing in the linearized error.

This observer is compatible with the bidomain model when it is rewritten in terms of $u^{(m)} = u^{(i)} - u^{(e)}$ and $u^{(e)}$ and when the gain correction is added to the equation on $u^{(m)}$. Always is [4], numerical illustrations including parameters estimation were performed to validate the strategy.

**Mumford-Shah functional** As discussed in [4], the observer design originates from an analogy with the dynamical behavior of contours tracking an object in a Mumford-Shah-based segmentation. The observer term corresponds to the shape derivative of the data-fitting term of the functional given in [CV91]

$$\mathcal{J}_u(u) = \int_{u>c_{th}} (y - C_{max}(y, u))^2 dS + \int_{u<c_{th}} (y - C_{min}(y, u))^2 dS,$$

that could also be seen as a functional of $\Omega_u = \{x \in \Omega, u > c_{th}\}$

$$\mathcal{J}_u(\Omega_u) = \int_{\Omega_u} (y - C_{max}(y, \Omega_u))^2 dS + \int_{\Omega \setminus \Omega_u} (y - C_{min}(y, \Omega_u))^2 dS,$$

as $C_{in}$ and $C_{out}$ could be rewritten easily in terms of $\Omega_u$.
1.4.2 A new state observer based on the topological gradient

Our objective in this work was to extend the state observer introduced in [4]. Indeed, this observer is limited in the sense that it can only correct the front that is already present in the simulation. This means that, if we do not have a propagation front to be associated with the data, we are not able to compute any observer correction. Our goal here is to circumvent this limitation, in order to track complex types of propagation patterns such as cardiac fibrillations corresponding to irregular and often very rapid heart rhythms. This work was published in the proceeding [5].

In order to improve the observer formulation, we intend to follow a strategy of increasing importance in shape optimization [BHR04] or level-set-based image segmentation [HKO07, HL09]. The central idea is to complement the required shape derivatives, used to modify the shape contours, by a topological derivative that represents the sensitivity of $J_u$ when removing a small part of the domain. Following [BHR04], we formally define the topological derivative of our functional by

$$d\mathcal{J}_u(\Omega_u) = \lim_{x \to 0} \frac{\mathcal{J}_u(\Omega_u \setminus B_{\rho,x}) - \mathcal{J}_u(\Omega_u)}{|\Omega_u \cap B_{\rho,x}|},$$

where $B_{\rho,x}$ denotes a ball of radius $\rho$ and center $x$. This gives in our case as in [HL09]

$$d\mathcal{J}_u(\Omega_u) = (y - C_{in}(y, \Omega_u))^2 - (y - C_{out}(y, \Omega_u))^2.$$

Then, we incorporate this topological derivative in our observer as a new reaction term, when we detect a topological difference between the target solution and the observer solution. In this respect, we follow a simple reasoning inspired from [HKO07]. Using illustrations, we show in [5] that we only want to act in the region where $d\mathcal{J}_u \times (u - c_{th})$ is positive. We infer that we should consider a topological gradient term of the form

$$\vec{\nabla}_{\text{top}} \mathcal{J}_u = (1 + \text{sign}(d\mathcal{J}_u \times (u - c_{th}))) d\mathcal{J}_u,$$

where the first term selects the region where we want to act while the topological derivative $d\mathcal{J}_u$ provides the direction (creation or destruction).

We illustrate this observer on a realistic configuration where we purport to track an atrial fibrillation starting from a physiologically healthy initial condition. In order to synthetically generate the atrial fibrillation, we use a standard S1-S2 protocol [CFPS14]. The location of the stimulation S2 is near the left pulmonary inferior vein. Indeed, the pulmonary veins are known to be prone to frequent re-entries. The S1 stimulus corresponds to a standard sinus stimulation (natural pacemaker of the heart) at t = 0 and t = 700 ms, corresponding to the cardiac cycle period considered here. The S2 stimulus is triggered at t = 290 ms. The results of the simulation of the target atrial fibrillation are displayed in Figure 1.6 (left column). The figure shows successive time steps between t = 300 ms (note that before the S2 stimulus the simulation is the same as in the healthy case) and t = 820 ms, and compares the reference solution with the results obtained with our observer, and with a direct model simulation undergoing S1 stimulation only like the observer. We see that the observer very quickly and accurately captures the complex nature of the fibrillation solution.

1.4.3 Application to real data

This part is in collaboration with the Inria team Carmen. This team is a partnership between Inria and IHU Liryc and is specialized in modeling of cardiac electrophysiology. IHU Liryc is a unique center in
<figure><img src="image_17.png" /><figcaption>Figure 1.6 – Atrial fibrillation: target solution (left), observer solution (middle) and direct simulation without data (right)</figcaption></figure>

the world which aims to understand and treat heart rhythm disorders better. This collaboration gave us access to a depolarization map of one patient extracted from Carto® system.

**Strategy sensitivity** Before applying it to real data, we evaluate the robustness of our data assimilation method when confronted with model simplifications. The model sensitivity is evaluated on the basis of synthetic data generated with a validated realistic atrial model that accounts for the specific atrial structure and fibers at endocardium and epicardium. The problem inverse is solved with simpler modeling components and this numerical analysis shows that, among other things, conductance anisotropy and knowledge of fiber distribution are essential to produce a model that matches observations, see [12]. Always in this article, we also work with real data but we do not present this result here. We will focus on another application detailed in the thesis of Antoine Gérard [Gér19] (in French).

**Data presentation and application** The goal of electro-anatomical mapping systems for clinicians is to give an anatomical representation of the chamber considered in the acquisition and to give electrical signal amplitude, signal complexity, local activation times... The functioning of the Carto® system used to register the data and the pre-processing steps to obtain data that can be used in practice are presented in the thesis of Antoine Gérard [Gér19] (in French). The data concerns the left atrium. The geometry and the depolarization map are presented in Figure 1.7-Top. One can see that even if the geometry is complete, the depolarization maps are partial. To deal with this problem, we consider instead of a constant gain parameter $\gamma$ a spatial gain $\gamma: x \mapsto 1_{\omega^{\text{obs}}}(x)$ where $\omega^{\text{obs}}$ is the domain in which the data are available allowing to act only if the data are available.

A direct application of our observer is its capability to smooth and complete the activation maps using
both corrections: the one based on shape gradient allow to smooth and to complete and the ones based on topological gradients to deal with the appearance of new fronts. Figure 1.7-Bottom illustrates the potential of the approach.

<figure><img src="image_18.png" /><figcaption>Figure 1.7 – Two views (left: leftanterior, right: posteroanterior) of depolarization maps of a patient left atrium. Top: Data extracted with Carto® system. Bottom: Map completion using sequential observers.</figcaption></figure>

**Remark 2.** The observer presented in Section 1.4.1 has also led to other applications. In the context of wild-fire spread, a version compatible with an eikonal equation was proposed and illustrated with synthetic data in [9] and then illustrated with real data in [14]. In this context, a version compatible with a Lagrangian front-tracking model was also developed [13]. In Chapter 2 (see Section 2.3), we also present a work directly inspired from the design of this observer in the context of tumor growth.

## Conclusion

This chapter presents my contributions to cardiac electrophysiology modeling. As for the modeling aspects, the origin of the very classical bidomain model was studied. A complete and exhaustive mathematical analysis of the periodic homogenization procedure leading to the macroscopic bidomain model when the mechanical deformations of the heart are neglected has been presented. An extension of the microscopic bidomain model by a microscopic bidomain model taking into account the deformations of the myocardium was also proposed and derived thanks to a formal homogenization inspired by our previous mathematical analysis. Very realistic simulations of cardiac electrophysiology validate our approaches. As for the aspects of the inverse problem, whose aim is to illustrate our desire to get as close as possible to clinical applications, the strategy – we developed and validated from a mathematical point of view – was applied to the data of a patient. This underlines that our approach has great potential. In my research project given in Chapter 4, I will present my future work in this area.
<figure><img src="image_19.png" /></figure>
# Chapter 2

## Contributions to tumor growth modeling

*The innovative work presented in this chapter has been published in four accepted articles [20, 19, 17, 24].*

### Introduction

The purpose of this chapter is to present my mathematical contributions to tumor growth modeling. This work has started since I arrived in Bordeaux in 2015. Developing models able to accurately describe tumor growth may help monitoring the disease evolution or even predicting the efficacy of different therapeutic strategies. Mathematical modeling of tumor growth has recently drawn interest from the medical community. Indeed, since the pioneering work on glioblastomas [SBMAJ03], the development of such reliable models has helped better understand biological mechanisms and clinical behaviors of several types of tumor such as lung metastases [CCJ+15], renal cell carcinoma [BBB+15], and glioblastoma [SLCFS14, SBMAJ03] or the effects of treatments [BAB+12, RKP+12, RSC+06]. Such applications are made possible mainly thanks to the development of realistic mathematical models which are personalized using patients' imaging.

My mathematical contributions to this application field can be summarized in the following points:
(1) **Select mechanistic models for tumor growth.** Unlike cardiac electrophysiology, where there is a reference model in the form of the bidomain model, there is no such consensus when modeling tumor growth. In my opinion, this can be explained by the fact that the unknowns of the system and the corresponding conservation laws and constitutive equations are still open questions due to the many uncertainties of how tumor growth works. Even though the roles of oxygen, nutrients, and angiogenesis (development of new vessels by the tumor itself) are becoming better understood, it can be difficult to translate them into mathematical equations, especially if one wants to obtain a parameterizable model with clinical data. In the Inria team MONC, we therefore decided to consider advection models for the temporal evolution of the densities of healthy and tumor cells.

(2) **Develop an efficient estimation strategy to adapt a tumor growth model to patients with meningioma.** Tumor growth models must be adapted to each individual case to generate predictive simulations for a given patient. In this context, abundant medical data – especially lesion volume, texture, and shape – can be used to identify model parameters. We develop an innovative strategy that combines a 0D reduction and population estimation approach to generate reliable predictions of the volume and the shape of tumors within a reasonable time frame. This strategy is illustrated using a large cohort of meningiomas (intracranial tumors). The results have been published in [17, 24].
(3) Adapt this strategy when the tumor shapes are time-evolving. The strategy described above is not efficient when the tumors have a time-evolving shape. To deal with this situation, we develop a sequential strategy. The main contribution is the development and the mathematical analysis of a Luenberger observer. This work has been published in [20].

These contributions are presented into three sections. Section 2.1 quickly presents the type of models used in the Inria team MONC and focus on a free growth model with only one density of tumor cells. The strategy developed to estimate its parameters very efficiently and using very few medical imaging is also presented. It corresponds to a part of the work published in [17]. In Section 2.2, the strategy is illustrated with a large cohort of intracranial meningiomas. Finally, in Section 2.3, we present the Luenberger observer to deal with time-evolving shapes whose the design is very inspired from the observers presented in Section 1.4.1 of Chapter 1.

## 2.1 Advection type models

### 2.1.1 General presentation

Let $\mathcal{B}$ be a bounded and smooth domain of $\mathbb{R}^3$. The tumor growth can be represented by mechanistic descriptions of the healthy and tumor cell densities evolution over time. We consider one (or more) tumor densitie(s) $P_i, \forall i \in [1; N_{TD}]$ – corresponding to multiple cell lines or at different status of the cell (proliferative, quiescent, necrotic) – and one healthy density $S$ such that $\sum_{i=1}^{N_{TD}} P_i + S = 1$. Historically in the MONC [RSC$^+$06, BRS$^+$09, BCG$^+$10] team, the models used are advection type models

$$ egin{cases} 
\partial_t P_i + \nabla \cdot (\vec{v} P_i) = \sum_{j=1}^{N_{TD}} f_{ij}(P_i, P_j), & \mathcal{B}, \quad \forall i \in [1; N_{TD}], \\ 
\partial_t S + \nabla \cdot (\vec{v} S) = 0, & \mathcal{B}. 
\end{cases} 
\quad (2.1) $$

One can see that we assume that $S$ is only advected at the velocity $\vec{v}$ due to cellular division of cancer cells. This velocity corresponds to the motion of tumor cells. This is clearly not true as there is a biological action of the tumor on the extracellular matrix but in our setting, this assumption means that the global evolution of the tumor is not modified by the extra-tumoral medium. The saturation assumption leads to

$$ \nabla \cdot \vec{v} = \sum_{i=1}^{N_T} \sum_{j=1}^{N_{TD}} f_{ij}(P_i, P_j). $$

To close the system, the Darcy's law is often considered $\vec{v} = -\nabla\pi$, where $\pi$ is the pressure that tumor cells exert on healthy cells. The functions $f_{ij}$ model the interactions between the different cells or the treatments impact. These complex models must be adapted to each individual case in order to produce predictive simulations for a given patient. In this context, the abundant available medical data – especially volumes, textures and shapes of the lesions – can be used to identify the model parameters.
### 2.1.2 Free growth tumor model

#### 2.1.2.1 Modeling

In this section, we present a growth model with only one tumor density (i.e. $N_{TD} = 1$). We assume that it is free growth meaning that there is no treatment. This implies that the function $f_{11}$ corresponds to the multiplication between $P$ and a growth rate denoted in that follows by $\tau_G$. The system reads

$$ \left\{ \begin{array}{rcl} \partial_t P + \nabla \cdot (\vec{v} P) & = & \tau_G(t, P)P, \quad \mathcal{B}, \\ \partial_t S + \nabla \cdot (\vec{v} S) & = & 0, \quad \mathcal{B}, \\ \nabla \cdot \vec{v} & = & \tau_G(t, P)P, \quad \mathcal{B}, \\ \vec{v} & = & -\nabla \pi, \quad \mathcal{B}. \end{array} \right. \qquad (2.2) $$

Different dynamics could be considered for the growth rate:

(1) a constant growth rate: $\tau_G(t, P) = a$,

(2) a time dependent decreased growth rate to model its consumption by the tumor: $\tau_G(t, P) = ae^{-bt}$, where $b^{-1}$ corresponds to the characteristic time at which the tumor growth capacity decreases and $a$ to the initial growth rate,

(3) or a spatial and time dependent growth rate: $\tau_G(t, P) = ae^{-b \int_0^t P(s, \cdot) ds}$ meaning that the consumption of the rate takes place only inside the tumor.

The two parameters $a$ and $b$ are assumed to be strictly positive constants. The modeling choices proposed here do not allow to obtain time-evolving shapes. To do so, one strategy consists in considering an initial growth rate distributed spatially meaning that $a$ is not anymore a constant. We will come back to this case in Section 2.3. It is clear than one could probably consider models for tumor growth closer to biology but the major advantage of these three models is that they depend only on one or two parameters: $a$ and $b$. We choose to keep the model simple to ensure that it can effectively be personalized for each patient using the sparse observations that are available in the clinical routine.

Concerning the boundary conditions (BC), it is quite easy to consider rigid constraints at the skull for example by considering that $\partial \mathcal{B} = \Gamma_c \cup \Gamma_f$, where $\Gamma_c$ is the boundary corresponding to the rigid constraint and $\Gamma_f$ the rest of the boundary. We then consider the following BC for the pressure

$$ \pi = 0 \text{ on } \Gamma_f \text{ and } \nabla \pi \cdot \vec{n} = \vec{v} \cdot \vec{n} = 0 \text{ on } \Gamma_c. \qquad (2.3) $$

We assume that the domain $\mathcal{B}$ is taken large enough such that the Dirichlet BC on $\Gamma_f$ has no significant influence on the growth.

Developing first equation on (2.2) yields

$$ \partial_t P + \vec{v} \cdot \nabla P = \tau_G(t, P)P(1-P). \qquad (2.4) $$

Let the initial condition of $P$ equal to $P(0, x) = 1$ inside the tumor and $0$ elsewhere, one can easily see that $P(t, X) = 1$ inside the tumor and $0$ elsewhere at any time. We define by $\Omega(t) = \{x, P(t, x) > 0.5\}$ the domain corresponding to the tumor at time $t$.

These models are simplified versions of the general model presented in Section 2.3.1 whose the existence and uniqueness is given in Propositions 2 and 3 under the hypothesis of Dirichlet boundary conditions for $\pi$ (see [20] for the proof). Simulation of these PDE systems could be done using a following split-step scheme described in [17].
### 2.1.2.2 ODE systems embedded in PDE systems

Let us denote by $V(t)$ the volume of the lesion $\Omega(t)$ at time $t$. We have

$$V(t) = \int_{\Omega} P(t, X) dX = \int_{\Omega(t)} dX$$

by definition of $P$. One can prove the following proposition using Reynolds and Green theorems, see [17].

**Proposition 1.** We have

$$ 
\begin{aligned}
V(t) &= V_0 e^{at}, \quad \text{if } \tau_G(t, P) = a, \\
V(t) &= V_0 e^{\frac{a}{b}(1-e^{-bt})}, \quad \text{if } \tau_G(t, P) = ae^{-bt}, \\
V(t) &= V_0 \left( 1 + a \frac{e^{(a-b)t} - 1}{a-b} \right), \quad \text{if } \tau_G(t, P) = ae^{-b} \int_0^t P(s, \cdot) ds \text{ and } a \neq b.
\end{aligned}
$$

The first case corresponds to an exponential increase of the volume. One can recognize in the second case the formula of the Gompertz model, a very classical growth model, see [BLB+14] for a full study and [Gom25, Cas34, Lai64, Lai65, Nor88] for validations of Gompertz model for tumor growth.

### 2.1.2.3 Estimation strategy

**Descriptive power** We will consider a cohort of $N_P$ patients with for each individual $j \in [1; N_P]$, $N_{obs}^j$ MRI or CT scan examinations acquired at different times $(t_k^j)_{k \in [1; N_{obs}^j]}$ are available. From these data, one can then extract – using for example a semi-automatic segmentation algorithm – the tumor volumes denoted $(\bar{V}_k^j)_{k \in [1; N_{obs}^j]}$, $\forall j \in [1; N_P]$. As PDE and ODE systems share the same parameters, they can be both estimated for a given patient $j$ by fitting the ODE systems to the tumor volumes denoted $(\bar{V}_k^j)_{k \in [1; N_{obs}^j]}$. This allows us to be able to deal with large cohort of real data.

Concerning theoretical identifiability, one can easily show that at least two (resp. three) measured volumes at two (resp. three) different times are necessary for the first (resp. second and third) case(s). However, only two or three medical imaging lead to overfitting and do not allow to integrate the fact that the data are noisy. One way to compensate for sparse sampling times and measurement uncertainties, is to consider a population approach allowing to constrain the variability of the parameters in the population of the $N_P$ patients. More precisely, one can for example consider mixed-effects approach, see Section B.2 of the introduction chapter. The initial volume instead of be fixed with the volume computed using the first MRI examination is then considered as a covariate of $\bar{V}_0^j$

$$V_0^j = \bar{V}_0^j(1 + e_0^j), \quad \forall j \in [1; N_P],$$

where $e_0^j \sim \mathcal{N}(0, \sigma_{V_0}^2)$. The parameters $a$ and $b$ are decomposed as

$$\log(\{a^j, b^j\}) = \log(\{a^{\text{pop}}, b^{\text{pop}}\}) + \{\zeta_a^j, \zeta_b^j\}, \forall j \in [1; N_P]$$

where $\{a^{\text{pop}}, b^{\text{pop}}\}$ are common for all the subjects and corresponds to fixed effects and $\{\zeta_a^j, \zeta_b^j\}$ – assumed to be centered – and are random effects. The Log-normal distribution allows to constrain the positivity
of the parameters. Concerning the measurement error, we consider it proportional

$$ \bar{V}_k^j = V^j(t_k)(1+e), \forall j \in [1; N_P]. $$

We assume that $e$ follows a centered Gaussian law. The standard deviation $\sigma_e$ is estimated and is used to define the 95% confidence interval: $[V^j(t)(1 - 1.96\sigma_e); V^j(t)(1 + 1.96\sigma_e)]$.

**Predictive power** We also develop a strategy to make prediction with only two medical imaging. More precisely, to answer the following question: using only two MRI examinations given at time $t_0$ and $t_1$, can we predict the volume and the shape of the tumor, at a later time $t_2$? It is a very difficult question because using only two time points, it is impossible to determine the behavior of the volume evolution for not constant growth rates, namely whether the volume evolution is linear, strictly convex or strictly concave. To answer this difficult question, we propose the following strategy. Let a patient denoted by $j$ for which only two time data are available, we assume that a population of patients with at least 3 times is available and we will use this population to estimate the parameters of the patient $j$ using once again the mixed-effects approach. Our strategy can be seen as a way to compensate the lack of available times of the given patient $j$ by an available population of $N_P - 1$ patients.

The estimation strategies have been done using the stochastic approximation expectation maximization (SAEM) algorithm [KL05] implemented in Monolix [Ant19]. For all results, the convergence of the algorithm has been reached. In [24], a validation of the strategy has been done using synthetic data.

## 2.2 Illustration on a large cohort of meningiomas

### 2.2.1 Meningiomas

Arising from arachnoidal cap cells, they account for a third of primary adult intracranial tumors [BGA$^{+11}$, OGF$^{+13}$]. They become a major public health issue: (1) an annual increasing incidence ranging from 2-4.5% has been reported [BEB$^{+18}$, BGA$^{+11}$, OGF$^{+13}$] and (2) one third of these tumors are asymptomatic at the time of diagnosis [BEB$^{+18}$, BGA$^{+11}$, HHI$^{+09}$, NFN$^{+05}$, OGF$^{+13}$]. While most meningiomas are slow-growing tumors, some may exhibit important growth over a short period of time. Recent data have shown that asymptomatic meningiomas are observed in 1% to 2% of the population, especially in women over 60 years of age [JOV$^{+11}$, KNP$^{+04}$, VIT$^{+07}$]. Surgery or radiation therapy remain the cornerstone of the treatment for meningiomas, and, following the diagnosis of asymptomatic lesions, clinicians struggle with the decision of either immediately treating the lesion using these invasive procedures, or simply planning periodical radiological examinations [BV97, GXCI$^{+19}$, IMM$^{+19a}$, KKC$^{+18}$, MCZ$^{+17}$, RRBP18, ZLJ$^{+15}$]. In most cases, the frequency of such a monitoring is based on empirical strategies [Cou06, HHI$^{+09}$, IMM$^{+19b}$, JKK$^{+11}$, MCZ$^{+17}$, NRM$^{+03}$, NNF$^{+11}$, OKSL11, RRBP18, SRA$^{+10}$, ZAD$^{+08}$, ZLJ$^{+15}$]. This is not optimal due to the potential toxicity of repeated Gadolinium infusions [LDAW18] and from a medico-economic point of view. Consequently, predicting the progression of asymptomatic lesions is paramount to inform decision-making. However, to our knowledge, only one predictive tool [BBRRH$^{+19}$, LKP$^{+17}$, LPPK17] is available to assess the growth rate of untreated meningiomas and help clinicians in selecting which patients to treat or simply monitor. This tool was developed based on statistical predictions using various clinical and imaging data collected at the time of diagnosis. Our objective is to use mathematical modeling of the underlying physiopathological processes to offer an accurate alternative in terms of prediction of growth.
<figure><img src="image_20.png" /><figcaption>Figure 2.1 - Left: MRI T1 of 4 patients with tumor boundaries in red. Right: Anatomical view of a meningioma with skull, brain and meningioma geometries delineated in 3D from the MRI examination.</figcaption></figure>

### 2.2.2 Clinical data

**Study design and participants** - In this retrospective study, we analyze the clinical, radiological and volumetric data of all the consecutive patients who were referred to the neurosurgical consultation of H. Loiseau between December 2013 and August 2021 for an incidental meningioma. All the patients had regularly empirically defined clinical and radiological follow-up. Treatment (surgery, stereotactic or fractioned radiotherapy) was proposed according to volumetric progression and/or neurological symptoms as done in daily practice and was approved by the tumor board.

**Procedure and variable extraction** - All the imaging (MRI scans with 3D millimetric T1 WI sequence with gadolinium or CT-scans with contrast product injection if contraindications to MRI) were analyzed with the Sophia Radiomics® software which allows 3D volumetric measures by a semi-automatically segmentation algorithm. Illustrations of the segmentation process are given in Fig 2.1. Using this tool, we have extracted for each patient of the cohort the volumes and the shapes of the tumor at each time. For very few of them (9 patients), the skull has also been segmented to perform 3D simulations.

**Inclusion and exclusion criteria for mathematical modeling** - Only tumors with at least three imaging at different time points were included. Exclusion criteria were missing baseline radiological data, cyproterone, nomegestrol and chlormadinone acetate treatment at the time of meningioma diagnosis, type 2 neurofibromatosis, and a statistically significant decrease in volume over time. During the study period, 294 consecutive patients with $N_P = 333$ tumors met the inclusion criteria.

**Clinical questions** - Using these data, we want to answer three questions:

(1) Descriptive power: Can we model the dynamics of meningioma volume evolution?

(2) Classification: Can we classify the meningiomas depending of their agressiveness?

(3) Predictive power: Using only 2 MRI examinations, can we predict the tumor volume and the shape at a later time?
<table><caption>Table 2.1 – Descriptive power: Fit performances according to AIC and BIC criteria, MSE and errors of the four volume models.</caption><thead><tr><th>Model</th><th>AIC</th><th>BIC</th><th>1/N<sub>P</sub> &sum;<sub>j=1</sub><sup>N<sub>P</sub></sup> MSE<sup>j</sup></th><th>1/N<sub>P</sub> &sum;<sub>j=1</sub><sup>N<sub>P</sub></sup> E<sup>j</sup><sub>L<sub>1</sub></sub></th><th>1/N<sub>P</sub> &sum;<sub>j=1</sub><sup>N<sub>P</sub></sup> E<sup>j</sup><sub>L<sub>&infin;</sub></sub></th></tr></thead><tbody><tr><td>Linear</td><td>1279</td><td>820</td><td>0.085 [0.056;0.14]</td><td>0.066 [0.045;0.11]</td><td>0.14 [0.091;0.24]</td></tr><tr><td>Exponential</td><td>983</td><td>1002</td><td>0.075 [0.047;0.13]</td><td>0.060 [0.038;0.10]</td><td>0.13 [0.077;0.22]</td></tr><tr><td>Power</td><td>802</td><td>829</td><td>0.074 [0.045;0.12]</td><td>0.056 [0.035;0.10]</td><td>0.12 [0.076;0.18]</td></tr><tr><td>Gompertz</td><td>793</td><td>819</td><td>0.066 [0.044;0.11]</td><td>0.053 [0.036;0.089]</td><td>0.11 [0.072;0.18]</td></tr></tbody></table>

### 2.2.3 Descriptive power

Concerning the volume evolution, we will consider the Gompertz model ($\tau_G(t) = ae^{-bt}$) and the exponential model ($\tau_G = a$). For the results with the spatial growth rate on a small cohort (40 patients), we refer to [17]. To answer the **first question**, the two models are also compared to two other very classical growth models which are not related to a PDE system:

* Linear model: $V(t) = V_0 + at$ ($V_0 > 0, a > 0$),

* Power model: $V(t) = (V_0^{1-\gamma} + a(1-\gamma)t)^{\frac{1}{1-\gamma}}$, solution of $V' = bV^{\gamma}$ ($V_0 = V(0) > 0, a > 0, 0 < \gamma < 1$).

To estimate the parameters of the four 0D models, the population approach presented in Section 2.1.2.3 is applied. The goodness of fit are compared by the Akaike (AIC) and the Bayesian information criterion (BIC). For each meningioma $j$, we also compute the errors between estimated and real volumes calculated according the relative $L^1$ and $L^\infty$ norms and the Mean Squared Error (MSE):

$$E_{L^1}^j = \frac{1}{N_{\text{obs}}^j} \sum_{k=1}^{N_{\text{obs}}^j} \frac{|V^j(t_k) - \bar{V}_k^j|}{\bar{V}_k^j}, \quad E_{L^\infty}^j = \max_{k=1 \cdots N_{\text{obs}}^j} |V^j(t_k) - \bar{V}_k^j|, \quad \text{MSE} = \frac{1}{N_{\text{obs}}^j} \sum_{k=1}^{N_{\text{obs}}^j} \frac{(V^j(t_k) - \bar{V}_k^j)^2}{(\bar{V}_k^j)^2}.$$

Table 2.1 gives AIC, BIC and errors. The Power and Gompertz models gave similar results with fairly close AIC and BIC and no significant difference between errors calculated by the 3 metrics. The Power and Gompertz models gave significantly smaller errors. Because the Gompertz model gave the best results in all the metrics used, it was selected to model the meningiomas growth. This model has also the advantage to be compatible with one of the PDE system allowing to test the availability to describe the shape of the meningiomas. The standard deviation $\sigma_e$ of the proportional error was estimated at 14.6%. Figure 2.2 depicts six representative growth curves of meningiomas fitted with the four models. Note how the linear model (resp. the exponential model) does not fit the data for ID 120 (resp. for ID 305). ID 50 and ID 53 are two tumors with very close initial volumes (9.2 and 9.9 cm³ resp.) at first imaging and with the same length of follow-up but with different growth pattern over time and thus for ID 53, only Gompertz model fits well the data. For ID 50 and ID 164, the 4 models are similar. ID 12 had similar dynamics to ID 53 despite a very smaller initial volume and once again, only Gompertz model fits well the data.

### 2.2.4 Classification

Once the parameters estimated, the estimated volume $V$ is known for all times using explicit formulas of the models. We then estimate the relative growth-rate (RGR) during the $X$ months before and after as

$$\text{RGR}_{-X}^j(t) = \frac{V^j(t) - V^j(t-X)}{V^j(t)} \quad \text{and} \quad \text{RGR}_{+X}^j(t) = \frac{V^j(t+X) - V^j(t)}{V^j(t)}, \quad \forall t \ge 0.$$
<figure><img src="image_21.png" /><figcaption>Figure 2.2 – Descriptive power - Six illustrative examples of tumor growth modeling with the four models.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 2.3 – Classification - Hierarchical clustering results.</figcaption></figure>
<table><caption>Table 2.2 – Predictive power - Relative error $E_{L1}^{t2}$ and relative bias $B_{L1}^{t2}$ (see Eq. (2.5)).</caption><thead><tr><th>Error (mean of 33 patients)</th><th><i>E</i><sub>L<sub>1</sub></sub><sup>t<sub>2</sub></sup><br>9.9%</th><th><i>B</i><sub>L<sub>1</sub></sub><sup>t<sub>2</sub></sup><br>0.007</th></tr></thead><tbody><tr><td>Gompertz model</td><td>12.5%</td><td>0.024</td></tr><tr><td>Linear model</td><td></td><td></td></tr></tbody></table>

To estimate if at a given time $t$, the RGR was increasing or decreasing, the following growth rate ratio (GRR) was computed $GRR_X^j(t) = \frac{RGR_{+X}^j}{RGR_{-X}^j}$. If $GRR_X^j$ is strictly inferior to 1, it means that the growth of the meningioma slows down, if $GRR_X^j$ equals to 1, the growth is linear and if $GRR_X^j$ is strictly superior to 1, the growth of the meningioma accelerates. One can easily prove that when volume follows a Gompertz model, the relative growth rates and the growth ratio are decreasing functions over time. To answer the **second question**, a hierarchical clustering of the meningiomas was done using $RGR_{+6}$ and $RGR_6$. The value of 6 months has been selected based on two facts: (1) our ability to well estimate the quantities at 6 months has been validated on synthetic data (see [24]) and (2) the period of time is sufficiently long to detect a change in growth rate. The dendrogram obtained was computed using the Ward.D2 agglomeration method of the hclust function of R. The number of clusters is selected using the dendrogram. Figure 2.3 gives the results of the classification. At the first examen time, the hierarchical clustering identifies 4 clusters – that could name – as

* **Cluster 1** - slowing-growth pattern meningiomas: $GRR_6 = 0.3$ [0.2; 0.4],

* **Cluster 2** - linear pattern meningiomas: $GRR_6 = 1.0$ [0.9; 1.0],

* **Cluster 3** - accelerating pattern meningiomas: $GRR_6 = 1.1$ [1.1; 1.2] and $RGR_{+6} = 19.4\%$ [14.3%; 28.8%],

* **Cluster 4** - very accelerating pattern: $GRR_6 = 1.6$ [1.4; 1.9] and $RGR_{+6} = 79.1\%$ [53.7%; 102.8%]. and at the last examen time, it identifies 3 clusters – that once again could name –

* **Cluster 1** - slowing-growth pattern meningiomas: $GRR_6 = 0.3$ [0.2; 0.5],

* **Cluster 2** - linear pattern meningiomas: $GRR_6 = 0.9$ [0.9; 1.0],

* **Cluster 3** - and accelerating pattern meningiomas: $GRR_6 = 1.2$ [1.1; 1.3].

The clinical and radiological features of the meningiomas by cluster have been studied by our collaborating clinicians. We only give a small summary to show how our work could be use useful but we refer to [24] for more details. We observe that age at diagnosis and the initial volume of the tumors significantly increased from cluster 4 to cluster 1 (from 56.9 to 64.7 years-old and from 0.34 to 1.43 cm³ respectively). We also observe a significantly increasing proportion of patients who have had a cranial radiotherapy before the diagnosis of meningioma from the cluster 4 to 1 (from 26.3% to 3.4% respectively). We did not find any correlation between the proportion of meningiomas with a hypersig- nal FLAIR, calcifications or edema and the clustering. Consistently, we do not observe any correlation between the clustering and the AIMSS risk score. We observe a significantly increasing proportion of treated meningiomas from cluster 1 to 4 (from 3.4% to 68.4% respectively) without any significant difference between the delay to treatment and volume at treatment. Among the meningiomas operated on, we observe a significantly increasing proportion of WHO grade II tumors from cluster 1 to 4 (from 0% to 64% respectively). The main perspective of this research is to try to predict the cluster using only the initial imaging features through radiomics approach. This work is in progress in the context of the PhD of Virginie Montalibet.
<table><caption>Table 2.3 – Predictive power - Sørensen-Dice coefficient SDC₂ (see Eq. (2.6))</caption><thead><tr><th>Pat</th><th>2</th><th>4</th><th>7</th><th>8</th><th>9</th><th>11</th><th>18</th><th>20</th><th>22</th><th>23</th></tr></thead><tbody><tr><td>PDE model</td><td>95%</td><td>93%</td><td>83%</td><td>87%</td><td>65%</td><td>85%</td><td>86%</td><td>90%</td><td>85%</td><td>87%</td></tr><tr><td>Redist.</td><td>94%</td><td>93%</td><td>74%</td><td>69%</td><td>52%</td><td>82%</td><td>86%</td><td>87%</td><td>66%</td><td>88%</td></tr></tbody></table>

### 2.2.5 Predictive power

Concerning the predictive power, the population strategy presented at the end of Section 2.1.2.3 has been applied in [17] on 33 patients for the volume prediction and on 10 patients for the shape prediction. Extend this database is one of the perspective of the PhD of Virginie Montalibet. The two main issues concern the optimization of the PDE resolution and the development of an automatic strategy to segment the skull.

To validate the strategy, we first compute the relative error $E_{L^1}^{t_2}$ and the relative bias $B_{L^1}^{t_2}$ between the predicted volume at $t_2$ and the measured volume given by

$$E_{L^1}^{t_2} = 100 \sum_{j=1}^{N_P} \frac{|V^j(t_2) - \bar{V}_2^j|}{\bar{V}_2^j} \quad \text{and} \quad B_{L^1}^{t_2} = \sum_{j=1}^{N_P} \frac{V^j(t_2) - \bar{V}_2^j}{\bar{V}_2^j}. \qquad (2.5)$$

Regarding the results given in Table 2.2, the prediction strategy gives very interesting results. Indeed, the results are excellent: using the ODE system, we can predict the volume at $t_2$ with an error inferior to 10% i.e. with the same order of magnitude of the segmentation error. Concerning the bias, it is smaller using the ODE system than linear model.

Second, concerning the shape, we will compare then using the Sørensen-Dice coefficient. To do so, we denote by $\bar{\Omega}_2^j$ the segmented 3D domain using the MRI at time $t_2$ for a given patient j and we denote by $\Omega_2^j$ the 3D domain corresponding to $\Omega_2^j = \{x \in \mathcal{B}, P^j(t_2, x) > 0.5\}$, where $P^j$ is the solution of the PDE system (2.2) for a given patient j. The Sørensen-Dice coefficient is defined as

$$\mathrm{SDC}_2^j = 100 \frac{2 |\Omega_j^i \cap \bar{\Omega}_2^j|}{|\Omega_j^i| + |\bar{\Omega}_2^j|}. \qquad (2.6)$$

On can easily see that for an individual meningioma SDC₂^j is in [0, 100]% and that if SDC₂^j is close to 100% (resp. to 0%), the predicted shape is close (far away) to the measured shape at $t_2$. We compare the solutions of the PDE model with the shape obtained by doing an expansion of the shape of the initial tumor using a redistancing process. To do that, we consider the level set of the initial shape with a volume equals to the volume predicted using the linear 0D model. Table 2.3 gives the results. One can see that the shape of the tumor is very well estimated with a mean of 85% for PDE systems and 79% for the linear model. The results are very encouraging. Figure 2.4 shows the shape for 10 patients. The limitation of the linear model (which is not able to take into account the boundary of the skull) can be seen on the results of Patients 7, 8, 9 and 22. Concerning the PDE system, one can see that – expected Patient 9 – the results are very encouraging. Concerning Patient 9, its shape evolves in time. It is one of the main limitation of the PDE model – presented in Section 2.1.2.1 – to not be able to integrate this phenomenon. In the next section, we will explain the approach that could be applied in this context.
<figure><img src="image_23.png" /><figcaption>Figure 2.4 – Predictive power - 3D views of the solution of PDE system (red) compared to the MRI segmentation (blue) at the last available time for 10 patients. The supportive arachnoid tissue (skull) is represented by the grey surface.</figcaption></figure>
## 2.3 Tumor model with time-evolving shapes

In this section, we present work published in [20] that aims to address the estimation of tumor evolution when they have a nontrivial shape. This estimation provides crucial information, for example, in brain or lung tumors, since the tumor may affect vital functions or structures of the diseased organ. Our strategy is based on the development of a joint state and parameter observer as presented in Algorithm (6) of the introduction chapter. Section 2.3.1 presents the model – which is a generalization of System (2.2) – and the existence and uniqueness results. The state observer is then given and studied in Section 2.3.2. Results of existence and uniqueness and convergence theorems are also given. In Section 2.3.3, the strategy to use in practice is presented. Numerical illustrations are given.

### 2.3.1 Tumor model and analysis

We consider a generalization of System (2.2) (corresponding to $\tau_G^0(t, P) = \tau_G^1(t, P) = \tau_G(t, P)P$) and with Dirichlet boundary conditions for the pressure

$$ egin{cases} 
\partial_t P + \nabla \cdot (\vec{v} P) = \tau_G^0(t, P), & \mathcal{B}, \\ 
S + P = 1, & \mathcal{B}, \\ 
\nabla \cdot \vec{v} = \tau_G^1(t, P), & \mathcal{B}, \\ 
\vec{v} = -\nabla \pi, & \mathcal{B}. 
\end{cases} 
\quad (2.7) $$

The functional $\tau_G^0$ describes the tumor evolution and the functional $\tau_G^1$ the induced pressure due to tumor expansion. We make the following assumptions on these two functionals.

**Assumption 2** (Assumptions on the functional mappings $\tau_G^i, i = 0, 1$).

* For $i = \{0, 1\}$, $\tau_G^i$ are uniformly bounded and smooth functionals which map $\mathbb{R}^+ \times L^\infty(0, T; H^s(\mathcal{B}))$ into $L^\infty(0, T; H^s(\mathcal{B}))$ for any $s \ge d/2$.

* For any $P \in L^\infty(0, T; H^s(\mathcal{B}))$, the function $g_P$ is defined on $\mathbb{R}^+ \times \mathcal{B}$ by

$$ g_P(t, x) = \tau_G^0(t, P(t, x)) - P\tau_G^1(t, P(t, x)). $$

The function $g_P$ is such that for any $P \in L^\infty(0, T; H^s(\mathcal{B}))$, $g_P$ belongs to $L^\infty(0, T; H^s(\mathcal{B}))$. Moreover if $P(t, \cdot)$ compactly supported in $\mathcal{B}$, then $g_P(t, \cdot)$ is compactly supported in $\mathcal{B}$.

The following propositions study the existence and uniqueness in case of Dirichlet boundary conditions on $\pi$ ($\Gamma_c = \emptyset$ in Equation (2.3)).

**Proposition 2.** Assume that the initial condition of $P$ belongs to $H^s(\mathcal{B})$ such that its support is compactly embedded in $\mathcal{B}$. Then there exists $T$ small enough such that Problem (2.7) coupled with $\pi = 0$ on $\partial\mathcal{B}$ admits a unique solution in $L^\infty(0, T; H^s(\mathcal{B}))$ whose support is compactly embedded in $\mathcal{B}$.

**Proposition 3.** Let $s \ge 5/2$ and assume that the initial condition of $P$ belongs to $H^s(\mathcal{B})$ such that its support is compactly embedded in $\mathcal{B}$, and assume in addition that its values belong to the interval $[0, 1]$. Let $\mathcal{L}_{\tau_G}^0$ satisfy Assumption 2, and set

$$ \forall u \in L^\infty(0, T; H^s(\mathcal{B})), \quad \tau_G^0(t, u) = \tau_G^1(t, u) = u\mathcal{L}_{\tau_G}^0(t, u). $$
For $i = 0, 1$, it is obvious that $\tau_G^i$ satisfies Assumption 2, and let $T$ be such that the unique solution $P$ of (2.7) exists in $L^\infty(0, T; H^s(\mathcal{B}))$, and such that the support of $P$ is compactly embedded in $\mathcal{B}$ then for any $t \in (0, T)$, $0 \le P(t, x) \le 1$.

The proofs given in [20] rely on the hyperbolic properties of the tumor model, combined with Gagliardo-Nirenberg estimates in Sobolev spaces.

### 2.3.2 State observer: design and study

**Data** As previously, we consider that we have at our disposal MRI or CT scan examinations at different times. In the previous presented work, only the segmented volumes are used to estimate the parameters. In this work, we investigate how the shape of the tumor may be used to circumvent the uncertainties associated with the dynamical model. More precisely, we consider that the measurements (denoted by $y$) are time-continuous sequences of tumor masks, namely, maps taking essentially two different values inside and outside the tumor region, up to some perturbations and regularization across the front.

**Context** We consider the particular case: $\tau_G^0(\cdot, P) = \tau_G^1(\cdot, P) = \mathcal{L}_{\tau_G}^0(\cdot, P)P$ where $\tau_G^i, i = \{0, 1\}$ verify Assumption 2. The target system reads

$$ 
\left\{
\begin{aligned}
\partial_t P + \nabla \cdot (\vec{v} P) &= \mathcal{L}_{\tau_G}^0(t, P)P, & \mathcal{B}, \\
\nabla \cdot \vec{v} &= \mathcal{L}_{\tau_G}^0(t, P)P, & \mathcal{B}, \\
\vec{v} &= -\nabla \pi, & \mathcal{B}, \\
P(0, \cdot) &= P_{\diamond}^0 + \zeta_P, & \mathcal{B},
\end{aligned}
\right.
\qquad (2.8) $$

where $\zeta_P$ is the initial condition uncertainty that we would like to control with the observer system

$$ 
\left\{
\begin{aligned}
\partial_t \hat{P} + \nabla \cdot (\vec{v} \hat{P}) &= \mathcal{L}_{\tau_G}^0(t, \hat{P}) \hat{P} + \gamma g_P(y, \hat{P}), & \mathcal{B}, \\
\nabla \cdot \vec{v} &= \mathcal{L}_{\tau_G}^0(t, \hat{P}) \hat{P}, & \mathcal{B}, \\
\vec{v} &= -\nabla \pi, & \mathcal{B}, \\
\hat{P}(0, \cdot) &= \hat{P}_{\diamond}^0, & \mathcal{B}.
\end{aligned}
\right.
\qquad (2.9) $$

In this work, we consider that $g_P$ designs a consistent observer providing physically realistic solutions.

**Assumption 3.** We assume that $g_P$ satisfies the following assumptions:

* If $P$ is the solution of the target system (2.8), then $g_P(y, P) = 0$,

* For any $u \in L^\infty(0, T; H^s)$ for $s$ large enough,

$$ 
\left\{
\begin{array}{ll}
\text{If } u \le 0 \text{ then } g_P(y, u) \ge 0, & \\
\text{If } u \ge 1 \text{ then } g_P(y, u) \le 0. &
\end{array}
\right.
\qquad (2.10) $$

**Proposition 4.** Assume that $P_{\diamond}^0$ and $\zeta_P$ are such that

$$ 0 \le P_{\diamond}^0 \le 1, \quad 0 \le P_{\diamond}^0 + \zeta_P \le 1, \quad \text{in } \mathcal{B}. $$

Assume in addition that $\mathcal{L}_{\tau_G}^0$ is such that the solution to Problem (2.8) exists in $L^\infty(0, +\infty; H^s)$. Finally assume that $g_P$ satisfies Hypothesis 3 and is such that for any $\gamma > 0$, the solution to Problem (2.9) exists in
$L^\infty(0, +\infty; H^s)$. Then

$$0 \le \hat{P} \le 1, \quad \text{a.e. in } (0, +\infty) \times \mathcal{B}.$$

The proof is obvious thanks to Assumption (2.10).

In order to investigate the convergence of the state observer towards the target solution, we consider the linearization of the observer term $g_P$

$$g_P(y, \hat{P}) \approx -dg_P(y, P)(P - \hat{P}).$$

This approximation is valid when the difference between $\hat{P}$ and $P$ is assumed to be small (i.e. the uncertainty of the initial condition is assumed to be reasonable). We denote by $\hat{P}$ the solution to the following linearized state observer model in the sense of distribution:

$$\left\{
\begin{array}{l}
\partial_t \hat{P} + \nabla \cdot (\vec{v} \hat{P}) = \mathcal{L}_{\tau_G}^0(t, \hat{P}) \hat{P} + \gamma dg_P(y, P)(\hat{P} - P), \quad \mathcal{B}, \\
\nabla \cdot \vec{v} = \mathcal{L}_{\tau_G}^0(t, \hat{P}) \hat{P}, \quad \mathcal{B}, \\
\vec{v} = -\nabla \hat{\pi}, \quad \mathcal{B}, \\
\hat{P}(0, \cdot) = \hat{P}_{\diamond}, \quad \mathcal{B}.
\end{array}
\right.
\qquad (2.11)$$

According to Proposition 2 there exists $T$ small enough such that the unique solution to Problem (2.8) exists and belongs to $L^\infty(0, T; H^s)$, $s > 5/2$. In addition if $dg_P(y, P)$ maps $L^\infty(0, T; H^s)$ into itself for $s > 5/2$, thus for any $\gamma > 0$, there exists a unique solution to Problem (2.11) in $L^\infty(0, T^\gamma; H^s)$, $s > 5/2$ with $T^\gamma < T$. However we are not interested in the condition to ensure global existence. We assume in the following proposition that $\mathcal{L}_{\tau_G}^0$ and $g_P$ are well-chosen to ensure a global existence of the solution to the two problems (2.8) and (2.11) and we focus on the convergence proof.

**Proposition 5.** Assume also that $P_{\diamond}^0$ and $\zeta_P$ are such that

$$0 \le P_{\diamond}^0 \le 1, \quad 0 \le P_{\diamond}^0 + \zeta_P \le 1, \quad \text{in } \mathcal{B}.$$

Assume in addition that $\mathcal{L}_{\tau_G}^0$ is such that the solution to Problem (2.8) exists in $L^\infty(0, +\infty; H^s)$ and that $g_P(y, P)$ is such that for any $\gamma > 0$, the solution to Problem (2.11) exists in $L^\infty(0, +\infty; H^s)$.

* If there exist $\alpha > 0$ and $\beta > 0$ such that

$$\int_{\mathcal{B}} dg_P(y, P)(u)udx \ge \alpha \|u\|_{L^2(\mathcal{B})}^2 - \beta, \quad \forall u \in L^\infty(0, +\infty; H^s), \qquad (2.12)$$

then there exists $c > 0$ such that for any $t > 0$, for any $\gamma > 0$,

$$\| (\hat{P} - P)(t) \|_{L^2(\mathcal{B})}^2 \le \| (\hat{P} - P)(0) \|_{L^2(\mathcal{B})}^2 e^{-(\gamma\alpha - c)t} + \frac{\gamma\beta(1 - e^{-(\gamma\alpha - c)t})}{\gamma\alpha - c}. \qquad (2.13)$$

In particular if $\beta = 0$, for any $\gamma > \frac{c}{\alpha}$, $t \mapsto \|(\hat{P} - P)(t)\|_{L^2(\mathcal{B})}$ decreases exponentially towards 0.

* If there exists $\alpha > 0$ and $p \ge 1$ such that

$$\int_{\mathcal{B}} dg_P(y, P)(u)udx \ge \alpha \|u\|_{L^p(\mathcal{B})}, \quad \forall u \in L^\infty(0, +\infty; H^s), \qquad (2.14)$$
there exists $\gamma_{0,p}$ such that for any $\gamma > \gamma_{0,p}$, $t \mapsto \|(\hat{P} - P)(t)\|_{L^2(B)}^2$ decreases linearly towards 0. In addition, if $p = 2$, then $t \mapsto \|(\hat{P} - P)(t)\|_{L^2(B)}^2$ decreases linearly towards 0.

The proof given in [20] consists in the study of the error system as explained in the first chapter. Inspired from the observer presented in Section 1.4.1 of Chapter 1, we consider

$$g_P(y, \hat{P}) = \delta_{\{x \in B, \hat{P}(t,x) = c_{th}\}}(t) \|\nabla \hat{P}\| ((y - C_{max}(y, \hat{P}))^2 - (y - C_{min}(y, \hat{P}))^2), \quad (2.15)$$

where $c_{th}$ is the threshold constant allowing to compare $\hat{P}$ to the segmented MRI domain. We refer to Section 1.4.1 of Chapter 1 for all the other notations. Always in [20], we show the following proposition.

**Proposition 6.** Assume that the data $y$ is sharp enough and satisfies the assumption given in [4]. Assume that there exists $\varepsilon > 0$ to be assumed small such that at any time $t$ the solution $\hat{P}$ of (2.11) with $g_P$ defined by Eq. (2.15) is such that $\tilde{P} = P - \hat{P}$ has a compact support included in a tubular neighborhood of $\Gamma_P(t) = \{x, P(t, x) = c_{th}\}$ of width $\varepsilon$. Then for any $\gamma > 2c\varepsilon/\alpha$, for any $t > 0$

$$\|\tilde{P}(t)\|_{L^2(B)}^2 \le \|\tilde{P}(0)\|_{L^2(B)}^2 e^{-(\gamma\alpha/\varepsilon - c)t} + 4\frac{\varepsilon^2\tilde{C}}{\alpha}.$$

The proof of the proposition is inspired from [4] and is based on the Fréchet derivative formalism. Regarding the choice of the gain parameter, the value of $\gamma$ through the proposition is probably not optimal and only roughly evaluable in practice. In [Kri18] (in French), a practical study to choose $\gamma$ is proposed.

**Remark 3.** Concerning the above proposition, the second hypothesis is similar to assume that the error is localized to the tumor front which is coherent in our context, where data consists mainly of tumor segmentation. The proposition also shows that the observer decreases the error exponentially fast from $\varepsilon$ to $\varepsilon^2$ in $L^2$ norm. Since the observer term only acts on the front, one cannot show the exponential decay over time.

### 2.3.3 And in practice?

To use the following state observer in practice, some points need to be lifted:

(1) **Selection of a model for $\tau_G^0 (= \tau_G^1)$**. In [20], we consider $\tau_G^0(t, P) = ae^{-b \int_0^t P(s, \cdot) ds}$ in the case where the initial growth rate $a$ is spatially distributed allowing to consider time-evolving shapes. When it is spatially distributed, it is denoted $a^s$.

(2) **Dealing with the very few available data.** In practical situations, data is not a time-continuous sequence of images, but a sequence of very few images acquired at distant times. In order to calculate the correction terms at a given time $t$, the idea is to interpolate data. The easiest way is to use a linear interpolation. However, when the tumor volume satisfies a solvable ODE, it is more accurate to use it to interpolate data. A redistanciation process has been taken in place, see [20] for all details.

(3) **Strong correlation between $a$ and $b$.** A correlation between the parameters $a$ and $b$ in a population of tumors has been observed frequently in Gompertz model. In [VRF+20], a reduced Gompertz model is even introduced: $b$ is estimated as explained in Section 2.1.2.2 and $k = a/b$ as a constant parameter in the population. In this context, we cannot consider this kind of strategy as we do not have a large database at our disposal. We only fix $b$ using the OD estimation obtained assuming $a$ constant. Only the initial distribution of the growth rate $a^s$ is estimated using the sequential strategy.
<figure><img src="image_24.png" /><figcaption>Figure 2.5 – Synthetic data obtained at time $t_0 = 0$, $t_1 = 15$ and $t_2 = 25$. This simulation is obtained with $a^s$ taking a value $a^{\text{in}} = 0.3$ inside the synthetic vessel (white front), and $a^{\text{out}} = 0.08$ outside the vessel.</figcaption></figure>

(4) **Dealing with the fact that the shape of $a^s$ is often not known.** Considering that $a^s$ is spatially distributed allows to integrate structures in the modeling. For example, $a^s$ can take different values in white and grey matters of the brain, because tumor propagation is faster in white matter than in grey matter [SAM00, HAS07]. We could consider also fibers and vessels surrounding the tumor that may also lead to shape changes. However, the locations of fibers and vessels are not visible using medical imaging and the segmentation of white and grey matter contains uncertainties and is computationally costly. Our objective is then to develop a strategy in which we have no prior on the structure of $a^s$ and to estimate it. We develop a strategy in which $a^s$ is estimated on a grid. Each case of the grid corresponds to a value that has to be estimated. The chosen grid is less refined that the computation grid, as Kalman filter is a costly algorithm.

(5) **Extend the value of $a^s$ outside the tumor to make prediction.** A smoothing and thresholding strategy has been developed and is presented in [Kri18] (in French).

In [20], numerical simulations with synthetic data allow to validate the Luenberger observer and the joint state and parameter observer on cases where the shape of $a^s$ is known with very few medical imaging. The last numerical simulation concerns the estimation in a case where the shape of $a^s$ is unknown. In order to simulate the growth of a tumor along a vessel, a target simulation is launched with $a^s$ taking a higher value inside a synthetic vessel. The simulation is shown on Figure 2.5, with the vessel front represented in white. Let us note that the tumor expands in the direction of the vessel as expected. The parameters used for this simulation are $(b, a^{\text{in}}, a^{\text{out}}) = (0.02, 0.3, 0.08)$, with $a^{\text{in}}$ and $a^{\text{out}}$ the values of $a^s$ inside and outside the vessel. Our objective is here to predict the tumor shape and to compare the real value of $a^s$ to the one estimated by our model. Only three times are used for the estimation process.

The volume calibration – done after assuming $a$ constant – gives $(b, a) = (0.029, 0.132)$. We will use the estimated value of $b$ for the PDE system and we will use the value of $a$ as an initial prior value for $a^s$. On the top of the first column of Figure 2.6, the first grid considered for $a^s$ is represented. On the bottom, the solutions of the target model in red and of the observer model in green are compared at time $t = 25$. We notice that the grid 5x5 contains three types of zones. Grey zones are those which are not reached by the tumor: parameters are not corrected, and stay at initial value. Red and blue zones are
<figure><img src="image_25.png" /><figcaption>Figure 2.6 – Initial shapes of $a^s$ (first line) and comparison between the solution of the target model and the solution of the observer model at time $t = 25$ (second line). First column: initial grid. Second and third columns: evolution of the grid using the previous results.</figcaption></figure>

those where the parameter increase or decrease respectively. In these simulations, we choose $\lambda = 0.001$ and an error of covariance of 0.08. It is interesting to notice that the red zones are those located inside the vessel. It shows that it is possible to extract information on the structure of $a^s$, even when we have no *a-priori* on it. In order to improve the prediction, the grid is refined like this: red zones are refined and blue and grey zones are unified to form a single zone. The aim is to refine interesting zones only. Unifying other zones aims at reducing the number of parameters and have reasonable computation time. We obtain successively the grids represented on the top of the second and third columns of Figure 2.6.

Solutions of the target and the observer models at time $t = 25$ on these grids are represented on the bottom of the second and third columns of Figure 2.6. We notice as before that red zones are those located in the synthetic vessel. We also notice that the most refined grid is not the one that gives the best result. Each zone being correlated to its neighbors, this method cannot give a precise vascularization map. However, the method is efficient to estimate a global behavior of the shape of $a^s$.

In this example, the strategy has been validated with synthetic data but in a clinical context. As these results are very encouraging, we also apply on a brain metastasis whose shape evolves in time (see Figure 2.7-Top). The results – summarized here – are detailed in [Kri18] (in French). The question of our collaborator was: can we determine the area with a large growth rate and which explains why the shape of the tumor evolves in time? We also add one: can we predict the volume and the shape at a later time using three MRI examinations? Figure 2.7-Middle/Bottom gives the results. The strategy allows us
DATA - 4 MRI examinations during one year

$t_0 = 0$

$t_1 = 65$

$t_2 = 275$

$t_3 = 324$

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

RESULTS - using the 3 first MRI examinations

<figure><img src="image_27.png" /><figcaption>Figure 2.7 – Top: MRI examinations of a brain metastasis of a patient followed by Ana Ortiz de Mendivil Arrate of HM Hospitales in Spain at 4 different times (given in days) and the 3D segmented tumors. Bottom: Tumor rate estimation. Volume and shape prediction.</figcaption></figure>
to well describe and even predict the volume of the tumor. Concerning the shape, one can first see that our strategy allows to determine the shape of the growth rate allowing to make prediction of the tumor. This example on real data validates the proof of context of our strategy.

## Conclusion

This chapter presents my contributions to tumor growth modeling. The tumor models that we clas- sically use in the MONC team to model tumor growth are presented. Different strategies to fit these models to each patient using available medical data are presented, as well as the results obtained in a large cohort of meningiomas.

In all examples presented here, no treatment was considered, but our approaches may be valid in this context as well. As an example, we present in [19] a model that integrates chemotherapy by con- sidering two different cell densities ($N_{TD} = 2$ in System (2.1)). To parameterize it, lesion textures were considered in addition to volume data. The strategy is illustrated on a small cohort of 17 patients with EGFR-mutated non-small cell lung cancer treated with tyrosine kinase inhibitors (TKI) who relapsed after an initial response to treatment (enlargement of the primary tumor, appearance of aggressive metastases, deterioration of the patient's general condition). We identified a biomarker calculated thanks to our mathematical model that allows us to determine whether the patient will have a survival under TKI based on only three initial CT scans.

Always in the context of cancer therapies, an increasing amount of my research focuses on electro- poration. Electroporation is a complex phenomenon that occurs when biological tissue is subjected to electrical pulses that allow either the direct killing of cells in the target tumor area or the introduction of molecules into living cells. The next chapter is devoted to the presentation of my work on this innovative therapy.
<figure><img src="image_28.png" /></figure>
# Chapter 3

# Contributions to electroporation modeling

*The innovative work presented in this chapter has been published in six articles [15, 18, 21, 25, 29, 28].*

## Introduction

This chapter is dedicated to the innovative medical therapy called electroporation. Electroporation consists in applying pulsed electric fields of high intensity and short duration (usually several hundred volts per centimeter for a hundred microseconds) to cells using catheters containing active electrodes (to apply the electric field) and often passive/inactive electrodes (which can be used to measure). The objective to create defects – called pores – at the surface of the cell membrane. Over the past two decades, this phenomenon has been studied in detail both *in vivo* and *in vitro*. Interestingly, this method has been shown to selectively affect cell membrane structures and their properties [TGR05], leading to therapeutic innovations. In particular, reversible electroporation allows the introduction of non permeant molecules (cytotoxic drugs such as bleomycin, DNA plasmids) into living cells without directly killing them [Mir01, AM10, Geh08, MSB<sup>+12</sup>], whereas at higher electric pulse intensity or duration, cells in the target region die by a nonthermal mechanism called irreversible electroporation (IRE) [DMR05]. IRE is then used as a novel nonthermal ablation modality: pulsed electric field ablation (PFA). It appears that PFA represents an interesting alternative to the usual nonsurgical ablation techniques – radiofrequency electromagnetic field (RFA) or cryotherapy – to treat deep-seated, inoperable tumors [MLR05, CEH<sup>+13</sup>] and, more recently, to perform cardiac ablations [KKI<sup>+19b</sup>, COJ<sup>+20a</sup>, RDN<sup>+21</sup>].

However, one of the main limitations to the use of electroporation in routine clinical practice arises from the technical difficulties associated with this therapy, in particular the difficulty of (1) accurately understanding cell death and (2) determining the treated zone. Mathematical modeling of the electroporation process at the cellular level could provide information that contributes to the understanding of this phenomenon, while modeling of the electroporation process at the tissue level could evaluate the efficacy of the treatment thanks to a prediction of the electric field distribution based on a mathematical modeling.

My mathematical contributions to this field are summarized in the following points:

(1) **Use mathematical modeling to better understand how the membrane becomes permeable.** The goal of this first contribution is to develop mathematical approaches to track cell membrane dynamics. In a first work, we develop a robust calibration strategy for the 4-electrode experimental setup – developed by our collaborators Tomás García-Sánchez and Luis M Mir, among others – that allows us to measure the impedance of cell myotubes subjected to electrical pulses [GSAL<sup>+15</sup>]. We then design and parameterize
an equivalent circuit model that allows to quantify the short-term behavior of membrane conductance few seconds after the pulse. This first work gives us informations about the global behavior of the cell membrane and in parallel, we work on the local behavior of the membrane. To do so, we develop a PDE model of membrane electroporation that combines membrane water content and transmembrane voltage, and an efficient numerical strategy that allows us to obtain numerical simulations of pore formation on the membrane. This work corresponds to the two articles [29] [25].

(2) **Compare and improve electroporation models at the tissue scale.** The objective of the second contribution is to develop mathematical approaches to model the electroporation process at the tissue scale. We first compare the solutions of the phenomenological electroporation models available in the literature. We also examine their calibrations using bioelectric measurements on rabbit livers performed by Sel et al. [SCB+05]. In a second work, we focus on the boundary conditions of the – very conductive and thin – passive electrodes. It is a difficult question because, unlike the active electrodes on which either the potential or the current is imposed – corresponding to Dirichlet or Neumann boundary conditions – for the inactive electrodes the value of the equipotential is an unknown of the PDE. This leads to a non-local constraint of the flux along the conductor interface, the so-called floating potential problem. We develop a strategy to address this problem in the context of electroporation. This strategy is compatible with all of the models of the literature. This work has been published in these two articles [15, 18].

(3) **Investigate the effects of irreversible and reversible electroporation on tumor growth.** The aim of this third contribution is to study the growth specificities of cancer cells subjected to a pulsed electric field. To do so we work on data on 3D tumor spheroids obtained by our collaborators Jelena Kolosnjaj, Muriel Golzio and Marie-Pierre Rols (IPBS, Toulouse University). Our strategy is based on an advanced advection model (compatible with System 2.1 of Chapter 2) and well-fitted estimation strategies and has been published in [21].

These contributions are presented in three sections. Section 3.1 presents work at the cellular scale. Section 3.2 presents the tissue-scale work. Finally, in Section 3.3, the study of the effects of electroporation on tumor spheroids is presented. The population-based estimation strategy for PDE systems presented in [28], is also introduced as it was applied on the same data.

**Remark 4.** The work corresponding to the first two contributions focuses on application in oncology but one can see that most of it could be adapted for use in cardiology as well. This point is further develop in my research project (see Chapter 4).

## 3.1 Membrane dynamics

One of the most important properties of a biological cell membrane is that it is a semipermeable layer, meaning that it allows only certain molecules to pass through. For example, a cell can regulate the concentration of certain ions inside it by letting in more or fewer ions from its environment through its membrane. The membrane of a cell consists largely of a lipid layer, which is why it behaves like a dielectric under the influence of an electric field. In contrast, the interior of the cell usually behaves like an electrolyte. When an electric field is applied to a cell, charge accumulation occurs at the membrane, resulting in a transmembrane voltage. This phenomenon has already been studied in Chapter 1 for cardiac cells, which have the peculiarity of being excitable and has been modeled with the microscopic bidomain model coupled with an ionic model, see Section 1.1.
When this voltage exceeds a certain threshold, the permeability of the cell membrane increases dra-
matically, so that molecules that previously could not enter the cell can now do so. This phenomenon is
known as electroporation. Although the phenomenon was discovered in the late 1960s, the way in which
membranes become permeable through the action of a strong electric field is still not well understood.
Mir's group has shown that the cell membrane remains permeable for several minutes after microsecond
pulses, which has led to a well-designed clinical protocol for electrochemotherapy [MOBP91, MNS<sup>22</sup>],
and Pakhomov et al. have shown similar results for nanosecond pulses [PKW<sup>07</sup>]. Mir's group has also
demonstrated the dual phenomenon of electroporation: rapid pore formation that closes quickly and a
long-term change in membrane permeability confirmed by molecular dynamics simulations [BDS<sup>12</sup>].

It is worth noting that despite extensive experimental evidence of increases in membrane permeability
following electrical pulses, there is still a lack of knowledge about the electrical properties of the mem-
brane in the minutes following the pulses. In this section, we present mathematical work that attempts
to fill these gaps. In a first subsection, we study the global behavior of the membrane based on the time
evolution of its resistance using bioimpedance measurements. In a second subsection, we focus on the
local behavior of the membrane by designing a microscopic bidomain model that integrates the state of
membrane evolution through a continuous phase order function $\phi: x \in \Gamma \mapsto \phi(x) \in [0, 1]$, where $\phi = 0$
represents the ideal pure lipid phase and $\phi = 1$ represents the ideal pure water phase.

### 3.1.1 Bioimpedancemetry

Bioimpedance measurement seems to be the method of choice to track the evolution of the electrical
properties of the membrane in the minutes following the pulses, since a large number of measurements
can be made. However, it is confronted with the difficult problem of electrode-solute interactions that
affect the measurements and their interpretations. Namely, when an electrode is immersed in an ion-
containing solute, a boundary layer called the Gouy-Chapman-Stern layer forms around the electrode.
In this layer, the ion concentrations - and thus the conductivity of the solute - differ significantly from
the concentrations in the bulk (far from the electrode). In other words, the electrode measurements
are greatly disturbed by this layer, and the conductivity of the medium can hardly be accurately re-
covered [Lvo12]. Moreover, the thickness of the layer increases (nonlinearly) with the electric current
injected into the electrode [Sch66, ITC<sup>13</sup>]. Therefore, a standard impedance measurement setup leads
to noisy and unusable results in biological experiments. Schwan et al. proposed an experimental setup
with 4 electrodes to limit the size of the solute-electrode interface of the measurements [SF68]. Roughly
speaking, 2 electrodes deliver the current and 2 others connected to a high input impedance measure
the potential. However, calibration of the system is far from trivial. Bao et al. proposed a well-designed
setup and the corresponding numerical calibration method for impedance spectroscopy of erythrocytes
in [BDS93]. The calibration method is based on impedance measurements of 3 samples whose electri-
cal properties are known. The resolution of a nonlinear system then allows the determination of the
calibration functions. García-Sánchez et al. proposed a different 4-electrode experimental setup that re-
quires a special calibration [GSAL<sup>15</sup>] because the method of Bao et al. leads to an unstable calibration.
Normally, impedance measurements of cells in culture are performed using very dense cell suspensions
to increase the cell volume to buffer volume ratio. By using fixed cells in a setup that does not require
cell separation (and respects the cell geometry and membrane characteristics), impedance measurements
can be performed under more physiological conditions. García-Sánchez et al. recently proposed such a
well-designed setup. In particular, they proposed an experimental setup for fast electrical impedance
spectroscopy measurements (1 spectrum per millisecond) to study the evolution of the electrical proper-
ties of adherent cell myotubes within milliseconds after the electroporation pulse [GSAL<sup>15</sup>].