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

HAL
open science

# Modeling individual disease trajectory: application to the prediction of treatment effect

Pierre-Emmanuel Poulet

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

## To cite this version:

Pierre-Emmanuel Poulet. Modeling individual disease trajectory: application to the prediction of treatment effect. Machine Learning [stat.ML]. EDITE de Paris, 2023. English. NNT: tel-04345196

HAL Id: tel-04345196

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

Submitted on 14 Dec 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling individual disease trajectory: application to the prediction of treatment effect

Author:

Pierre-Emmanuel POULET

Aramis lab, INRIA Paris, Institut du Cerveau, Inserm, CNRS, Sorbonne Université

Supervisor:

Stanley DURRLEMAN

Directeur de recherche INRIA

PhD Thesis for the degree of Doctor in Applied Mathematics at EDITE Paris (ED 130)

Date of defense: 2nd of October 2023
Modeling individual disease trajectory: application to the prediction of treatment effect

Author:
Pierre-Emmanuel POULET (Aramis lab, INRIA Paris, Institut du Cerveau, Inserm, CNRS, Sor-bonne Université)

Supervisors:
Stanley DURRLEMAN (Directeur de recherche INRIA)
Abstract

This thesis aims at developing statistical models to describe the progression of neurodegenerative diseases from longitudinal data, i.e. repeated measurements in time for every subject. Our work builds upon the Disease Course Mapping, a spatio-temporal disease progression model. This mixed-effect model has a hierarchical structure allowing to account for the global dynamics at a population level while also accounting for individualized disease trajectories. The variability between subjects is decomposed between two effects: a temporal reparametrization to construct a common disease timeline, and a spatial effect to describe the heterogeneity in how the disease presents itself. However it suffers from several limitations in addressing the challenges of disease modelling. First this model is embedded in a Riemannian manifold, where the family of possible trajectories is restricted by the metric. Secondly it only considers continuous manifold-valued observations, while neurodegenerative cohorts are filled with discontinuous data. Thirdly it lacks a structure to identify and characterize disease subtypes. We will address all those needs.

In the first part, after introducing the Disease Course Mapping model, we propose a model that loosens the choice of possible trajectories by learning the metric of the Riemannian manifold.

In the second part we tackle discrete data and extend the model to binary and ordinal markers of disease progression. This allows us to predict symptoms and model cognitive scales. We apply it to Parkinson's disease and derive clinical interpretation from the results.

The third part is dedicated to heterogeneity. We extend the progression model to account for disease subtypes by adding a mixture layer on top of the statistical model. The mixture model is able to jointly uncover latent classes while learning the typical progression pattern of each class. We showcase the model's ability to extract distinct trajectories of disease evolution in the Alzheimer's disease case.

The fourth part focuses on the application of disease progression to treatment effect models. We define a framework to learn a therapeutic effect using only treated patients by leveraging a disease model's predictions, which act as a synthetic control arm. Application to the dopaminergic treatment in Parkinson's disease shows great promise. We then design a new model to account for disease-modifying treatment effects, i.e. long-term modifications of the disease progression trajectory.

Overall this thesis provides several tools to extend the Disease Course Mapping potential applications, easing its transfer towards having meaningful clinical impact. The methods developed here are aimed at improving our understanding of neurodegenerative diseases, by identifying subtypes, including markers of various nature and taking into account the effect of treatment over these markers.
<figure><img src="image_4.png" /></figure>
# Résumé

Cette thèse vise à développer des modèles statistiques pour décrire la progression des maladies neu-
rodégénératives à partir de données longitudinales, c'est-à-dire des mesures répétées dans le temps
pour chaque sujet. Notre travail s'appuie sur le modèle Disease Course Mapping, un modèle de
progression de la maladie spatio-temporel. Ce modèle à effets mixtes a une structure hiérarchique
permettant de rendre compte de la dynamique globale au niveau de la population tout en tenant
compte des trajectoires individuelles de la maladie. La variabilité entre les sujets est décomposée
en deux effets : une reparamétrisation temporelle pour construire une chronologie commune de la
maladie, et un effet spatial pour décrire l'hétérogénéité de la manifestation de la maladie. Cepen-
dant, il présente plusieurs limitations pour relever les défis de la modélisation de telles maladies.
Tout d'abord, ce modèle est défini dans une variété riemannienne, où la famille de trajectoires
possibles est restreinte par la métrique. Deuxièmement, il ne considère que des observations con-
tinues à valeurs dans la variété, alors que les cohortes neurodégénératives sont remplies de données
discontinues. Troisièmement, il manque une structure pour identifier et caractériser les sous-types
de maladies. Nous répondrons à tous ces besoins.

Dans la première partie, après avoir présenté le modèle Disease Course Mapping, nous proposons
un modèle qui assouplit le choix des trajectoires possibles en apprenant la métrique de la variété
riemannienne.

Dans la deuxième partie, nous abordons les données discrètes et étendons le modèle aux mar-
queurs binaires et ordinaux de progression de la maladie. Cela nous permet de prédire les symptômes
et de modéliser les échelles cognitives. Nous l'appliquons à la maladie de Parkinson et tirons une
interprétation clinique des résultats.

La troisième partie est consacrée à l'hétérogénéité. Nous étendons le modèle de progression
pour tenir compte des sous-types de maladies en ajoutant une couche de mélange au-dessus du
modèle statistique. Le modèle de mélange permet d'identifier conjointement les classes latentes
tout en apprenant le schéma de progression typique de chaque classe. Nous présentons la capacité
du modèle à extraire des trajectoires distinctes de l'évolution de la maladie dans le cas de la maladie
d'Alzheimer.

La quatrième partie se concentre sur l’application aux modèles d’effet de traitement. Nous définissons un cadre pour apprendre l’effet thérapeutique avec uniquement des patients traités en utilisant les prédictions d’un modèle de progression de la maladie, qui agissent comme un bras de contrôle synthétique. L’application au traitement dopaminergique dans la maladie de Parkinson illustre le potentiel de la méthode. Nous concevons ensuite un nouveau modèle pour tenir compte des effets de traitement modificateurs de la maladie, c’est-à-dire des modifications à long terme de la trajectoire de progression de la maladie.

Dans l'ensemble, cette thèse propose plusieurs outils pour étendre les applications potentielles du
Disease Course Mapping, facilitant son transfert vers un impact clinique. Les méthodes développées ici visent à améliorer notre compréhension des maladies neurodégénératives, en identifiant les sous-types, en incluant des marqueurs de nature variée et en tenant compte de l'effet du traitement sur ces marqueurs.
Remerciements

Tout d'abord je voudrais remercier Stanley de m'avoir encadré, d'abord pendant mon stage de fin d'études et ensuite pendant ma thèse. Tu m'as permis d'évoluer dans le monde de la recherche en stimulant mon questionnement, en me laissant explorer mes propres solutions et en me ramenant sur le droit chemin quand je partais trop loin. Ensuite je voudrais remercier les anciens du labo : Igor, Raphaël, Tiziana, Simona, Juliana, Etienne, Arnaud, Sophie, qui sont tous partis depuis (sauf Juli). C'est grâce à vous que je suis venu à l'Institut, et même si nous n'avons pas tous eu l'occasion de beaucoup nous voir à cause du Covid qui a marqué le début de ma thèse, vous avez tous participé à me faire aimer cet endroit qui est honnêtement le meilleur cadre de travail que j'ai pu voir et avoir. Ce cadre il doit beaucoup aux gens, à l'ambiance amicale et ouverte qui règne dans l'équipe. Je vous remercie donc sincèrement tous d'être là au quotidien : Némo, Juliette, Cécile, Ravi, Matthieu, Charley, Rémy, Tristan, Vito, Lisa, Camille, Elise, Arya, Sophie, Maëlys. Merci à tous les co-thésards avec qui nous avons évolué dans les méandres des publications, de l'administratif et surtout des afterworks. Merci à tous les ingénieurs qui nous rappellent au quotidien qu'un code de recherche ce n'est qu'un brouillon. Merci à Ninon de t'occuper de tous nos petits soucis informatiques.

Je remercie en particulier toute l'équipe Leaspy, avec qui les discussions sont toujours sources d'un éternel questionnement sur nos modèles. Je remercie mes collaborateurs, en particulier Bruno qui a suivi mon travail depuis deux ans. Ensuite je remercie les stagiaires que j'ai eu la chance d'encadrer, Samuel et Paul, qui sont maintenant en thèse. Vous avez apporté de la nouveauté dans mon regard sur mon travail en général et vos publications sont la preuve de la pertinence de réflexions.

Voilà j'en ai fini avec Aramis, du moins pour l'instant parce que l'aventure continue avec vous pendant encore un an.

Maintenant je me dois de remercier ma famille pour avoir soutenu mon parcours depuis le plus jeune âge. J'ai une pensée énorme pour mes grands-parents maternels qui se sont beaucoup investis dans ma réussite, et ce doctorat est une fierté de plus pour moi, une façon de leur dire merci.

Merci à tous mes amis d'école, les expats du Volley en famille élargie, dont beaucoup font une thèse et connaissent ce chemin si particulier. J'aurais cette petite fierté d'avoir été le premier du groupe à obtenir ma thèse.

Enfin merci à Nell, qui me supporte depuis maintenant plus de cinq ans. Je ne sais pas ce que j'aurais fait sans toi.
<figure><img src="image_5.png" /></figure>
# Contents

<table><tr><td>Abstract</td><td>iii</td></tr><tr><td>Résumé</td><td>v</td></tr><tr><td>Remerciements</td><td>vii</td></tr><tr><td>Introduction</td><td>1</td></tr><tr><td>Motivation</td><td>2</td></tr><tr><td>Disease progression modelling</td><td>7</td></tr><tr><td>Parkinson's disease</td><td>9</td></tr><tr><td>Thesis goals</td><td>11</td></tr><tr><td>Manuscript overview</td><td>12</td></tr><tr><td>Introduction en français</td><td>13</td></tr><tr><td>Motivation</td><td>13</td></tr><tr><td>Modèles de progression de la maladie</td><td>19</td></tr><tr><td>Maladie de Parkinson</td><td>21</td></tr><tr><td>Buts de la thèse</td><td>24</td></tr><tr><td>Plan du manuscrit</td><td>24</td></tr><tr><td>1 The Disease Course Mapping model: limits and first extension</td><td>27</td></tr><tr><td>1.1 The Disease Course Mapping model</td><td>29</td></tr><tr><td>1.1.1 Notations and mixed-effect modelling</td><td>29</td></tr><tr><td>1.1.2 Generic disease course mapping framework</td><td>32</td></tr><tr><td>1.1.3 Bayesian model</td><td>36</td></tr><tr><td>1.1.4 Limitations</td><td>38</td></tr><tr><td>1.2 Geodesic Bending</td><td>38</td></tr><tr><td>1.2.1 Introduction and related work</td><td>39</td></tr><tr><td>1.2.2 Method</td><td>40</td></tr><tr><td>1.2.3 Experiments</td><td>43</td></tr><tr><td>1.2.4 Discussion on Geodesic Bending</td><td>50</td></tr><tr><td>1.3 Conclusion</td><td>50</td></tr></table>
<table><tbody><tr><td><b>2 Discrete data</b></td><td><b>51</b></td></tr><tr><td>2.1 Context: discrete data, symptoms and disease scales</td><td>53</td></tr><tr><td>2.2 Binary model</td><td>56</td></tr><tr><td>2.2.1 Equation</td><td>56</td></tr><tr><td>2.2.2 Application to symptom's prediction</td><td>56</td></tr><tr><td>2.3 Ordinal model</td><td>61</td></tr><tr><td>2.3.1 Equation</td><td>61</td></tr><tr><td>2.3.2 Experiments: Simulation study</td><td>62</td></tr><tr><td>2.3.3 Application to Parkinson's disease data</td><td>66</td></tr><tr><td>2.3.4 Discussion</td><td>77</td></tr><tr><td>2.4 Parallel with Item Response Theory</td><td>77</td></tr><tr><td>2.4.1 Fisher information</td><td>79</td></tr><tr><td>2.5 Conclusion</td><td>84</td></tr><tr><td><b>3 Subtyping</b></td><td><b>85</b></td></tr><tr><td>3.1 Context</td><td>86</td></tr><tr><td>3.2 Mixture of DCM models</td><td>88</td></tr><tr><td>3.3 Mixture of individual parameters</td><td>89</td></tr><tr><td>3.4 Experiments</td><td>90</td></tr><tr><td>3.4.1 Simulated data: Univariate model</td><td>90</td></tr><tr><td>3.4.2 Applications to Alzheimer's disease data</td><td>91</td></tr><tr><td>3.4.3 Experiment with mixture of individual parameters</td><td>95</td></tr><tr><td>3.4.4 Combination with DCM models on discrete data</td><td>97</td></tr><tr><td>3.5 Discussion</td><td>99</td></tr><tr><td>3.6 Conclusion</td><td>99</td></tr><tr><td><b>4 Treatment models</b></td><td><b>101</b></td></tr><tr><td>4.1 Introduction to treatment modelling</td><td>103</td></tr><tr><td>4.1.1 Modelling treatment effect</td><td>103</td></tr><tr><td>4.1.2 Disease-modifying treatment effect</td><td>105</td></tr><tr><td>4.2 Additive treatment effect model</td><td>106</td></tr><tr><td>4.2.1 Latent disease progression model</td><td>109</td></tr><tr><td>4.2.2 Treatment effect model</td><td>109</td></tr><tr><td>4.2.3 Estimation method</td><td>112</td></tr><tr><td>4.2.4 Experiments</td><td>113</td></tr><tr><td>4.2.5 Discussion</td><td>128</td></tr><tr><td>4.3 Piecewise-geodesic model</td><td>129</td></tr><tr><td>4.3.1 Method</td><td>129</td></tr><tr><td>4.3.2 Experiments</td><td>132</td></tr><tr><td>4.3.3 Discussion</td><td>134</td></tr><tr><td>4.4 Conclusion</td><td>137</td></tr><tr><td><b>5 Estimation</b></td><td><b>139</b></td></tr><tr><td>5.1 The MCMC-SAEM algorithm</td><td>140</td></tr><tr><td>5.2 Variants of the MCMC SAEM</td><td>143</td></tr><tr><td>5.2.1 MMCMC SAEM for mixture models</td><td>143</td></tr><tr><td>5.2.2 Alternating maximization algorithm</td><td>146</td></tr></tbody></table>
CONTENTS

<table><tr><td>5.2.3 Gradient maximization step</td><td>148</td></tr><tr><td>5.2.4 Discussion</td><td>150</td></tr><tr><td>Conclusion and perspectives</td><td>151</td></tr><tr><td>Publications and scientific communication</td><td>155</td></tr><tr><td>Software development</td><td>157</td></tr><tr><td>Appendix</td><td>159</td></tr><tr><td>Riemannian geometry</td><td>159</td></tr><tr><td>Disease Course Mapping formulas</td><td>162</td></tr><tr><td>Introduction</td><td>162</td></tr><tr><td>Model description</td><td>163</td></tr><tr><td>Bayesian model</td><td>164</td></tr><tr><td>Likelihood</td><td>166</td></tr><tr><td>Gradients</td><td>168</td></tr><tr><td>Fisher information</td><td>170</td></tr><tr><td>Bibliography</td><td>185</td></tr></table>
<figure><img src="image_6.png" /></figure>
# Introduction

## Contents of the chapter

<table><tr><td>Motivation</td><td>2</td></tr><tr><td>Disease progression modelling</td><td>7</td></tr><tr><td>Parkinson's disease</td><td>9</td></tr><tr><td>Thesis goals</td><td>11</td></tr><tr><td>Manuscript overview</td><td>12</td></tr></table>
Motivation

In the wake of medical and sanitary progress over the last centuries, our societies have gained considerably in life expectancy. With the advent of an ever-aging population, new healthcare challenges have risen, namely cancer and neurodegenerative diseases. Focusing on the latter, this category of pathologies encompasses a large spectrum of diseases, being defined by a dysfunction in the peripheral or central nervous system. Currently neurodegenerative diseases are the second leading cause of deaths worldwide<fnref n="1" />. It is estimated that one in two women and one in three men will be affected by it during their life. Dementia affects 1 person in 14 after the age of 65 and 1 person in 6 after 80 years old, Alzheimer's disease (AD) being the most frequent cause of dementia with about 70% of all dementia cases. The second most common neurodegenerative disease is Parkinson's disease (PD), with more than 12 millions patients in Europe. These figures are bound to increase in the years to come, as demographic previsions forecast an increase of the 65 and older from 18% to 24% in 2060 (French population statistics<fnref n="2" />). However these neurodegenerative diseases raise several issues in their study.

**Understanding the disease:** Contrarily to infectious diseases, where the cause is easily identifiable and the pathogen mechanisms can be studied, neurodegenerative diseases remain poorly understood for the most part. There is no bacteria or virus to blame, no clear debut such as an infection, no cure nor vaccine. Their pathophysiology have been based on symptoms and clinical assessments of cognitive or motor abilities, which are the result from advanced neurodegeneration. However the process of neurodegeneration started much earlier, undetected before it is too late. For instance in AD or PD, the disease can last up to decades, and evidence shows that the biological mechanisms responsible for neuronal damage started more than twenty years prior to diagnosis<fnref n="3" />. This hidden evolution of neurodegeneration has been the main reason why medical research has been lagging behind in the understanding of the diseases. Diagnosed patients only account for the late part of the disease, when most of the damage has been done. The study of the natural disease history therefore is hampered by a lack of subjects with early disease development, in a state which is called prodromal. The recruitment of prodromal individuals is particularly difficult, due to another missing stone in our knowledge of neurodegenerative diseases: the causes remain mostly unknown. Risk factors have been identified<fnref n="4" />,<fnref n="5" />, including environmental factors or genetics. However they only account for a part of the diseases cases, and they do not automatically imply that the individual at risk will develop the pathology, which is the nuance between a risk factor and a cause. We are often faced with the large majority of patients being affected spontaneously, which is referred to as the idiopathic form of the disease.

The study of those diseases have led to hypothetical models trying to shape our understanding of their progression. In Alzheimer's disease the famous amyloid cascade mechanism was proposed

<footnote n="1">¹J. Dumurgier and C. Tzourio, "Epidemiology of neurological diseases in older adults," en, *Revue Neurologique*, vol. 176, no. 9, pp. 642-648, Nov. 2020.</footnote>

<footnote n="2">²*Projections de population à l'horizon 2060 - Insee Première - 1320.*</footnote>

<footnote n="3">³R. B. Postuma, A. E. Lang, J. F. Gagnon, et al., "How does parkinsonism start? Prodromal parkinsonism motor changes in idiopathic REM sleep behaviour disorder," *Brain*, vol. 135, no. 6, pp. 1860-1870, Jun. 2012.</footnote>

<footnote n="4">⁴G. Livingston, J. Huntley, A. Sommerlad, et al., "Dementia prevention, intervention, and care: 2020 report of the Lancet Commission," *English, The Lancet*, vol. 396, no. 10248, pp. 413-446, Aug. 2020.</footnote>

<footnote n="5">⁵T. Nedelec, B. Couvy-Duchesne, F. Monnet, et al., "Identifying health conditions associated with Alzheimer's disease up to 15 years before diagnosis: An agnostic study of French and British health records," *English, The Lancet Digital Health*, vol. 4, no. 3, e169-e178, Mar. 2022.</footnote>
in<fnref n="6" />. This model highlighted the role of the β-amyloid protein, whose progressive accumulation in the brain led to neuronal death, ultimately causing dementia. This hypothetical model has shaped the research in the field since then, to the point where the aggregating amyloid proteins are the primary target of drug trials. The study of this cascade over observable aspects of the disease has led to another widely known hypothetical model in AD. Introduced in<fnref n="7" />, the authors proposed that the evolution of disease markers, such as the amyloid burden, neuronal death markers and cognitive decline, would follow a logistic progression from normal to abnormal. Similar models have also been proposed in PD, namely the Braak stages<fnref n="8" /> and more recently the Brain-first and Body-first hypothesis<fnref n="9" />. The crucial part in those descriptive models is to provide an *ordering* of markers degradation, which is then used to define stages in the disease. Staging is essential as we will see when discussing the development of a treatment.

**Defining the disease:** This fledgling understanding of the biological processes underlying the observable pathophysiology of the disease questions even the very definition of these diseases. Until quite recently the notion of disease was related to the clinical manifestations, i.e. symptoms, and the diagnosis was assessed by a neurologist. However the patient probably *had* the disease before the symptoms' onset. Therefore recent developments push for a biological definition of the diseases, based on objective measurements of biomarkers which would allow for a diagnosis prior to the clinical one. In AD such a definition relies on the amyloid presence in the brain, or the tau burden -another protein whose misfolding induces tangles and fibrils, toxic for the neurons-, and patients are thus described under an A/T/N (for Amyloid/Tau/Neurodegeneration) classification<fnref n="10" />.

**Embrace heterogeneity:** This last example of classification in AD patients echoes the wide range of possible disease manifestations. Oftentimes a neurodegenerative disease encompasses various subtypes, sometimes due to specific genetic variants. This heterogeneity in the disease spectrum also explains why it is hard to define the disease properly. Heterogeneity is unavoidable in some cases because studying smaller subgroups of patients would leave too few data to study the disease, and too few potential subjects to be recruited in a clinical trial. Therefore one always has to account for the shared part of the disease physiology as well as the differences between the subtypes. The distinction between subtypes is also not always straightforward, as there is not always a genetic variant or a separating marker. Many attempts have been made for disease subtyping in PD, for instance<fnref n="11" />,<fnref n="12" />,<fnref n="13" />, but lines between the subgroups are blurred and not always stable as the diseases

<footnote n="6">⁶J. A. Hardy and G. A. Higgins, "Alzheimer's disease: The amyloid cascade hypothesis," eng, *Science* (New York, N.Y.), vol. 256, no. 5054, pp. 184-185, Apr. 1992.</footnote>

<footnote n="7">⁷C. R. Jack, D. S. Knopman, W. J. Jagust, et al., "Hypothetical model of dynamic biomarkers of the Alzheimer's pathological cascade," English, *The Lancet Neurology*, vol. 9, no. 1, pp. 119-128, Jan. 2010.</footnote>

<footnote n="8">⁸H. Braak, K. D. Tredici, U. Rüb, et al., "Staging of brain pathology related to sporadic Parkinson's disease," en, *Neurobiology of Aging*, vol. 24, no. 2, pp. 197-211, Mar. 2003.</footnote>

<footnote n="9">⁹J. Horsager, K. B. Andersen, K. Knudsen, et al., "Brain-first versus body-first Parkinson's disease: A multi-modal imaging case-control study," *Brain*, vol. 143, no. 10, pp. 3077-3088, Oct. 2020.</footnote>

<footnote n="10">¹⁰C. R. Jack, D. A. Bennett, K. Blennow, et al., "A/T/N: An unbiased descriptive classification scheme for Alzheimer disease biomarkers," eng, *Neurology*, vol. 87, no. 5, pp. 539-547, Aug. 2016.</footnote>

<footnote n="11">¹¹R. S. Eisinger, C. W. Hess, D. Martinez-Ramirez, et al., "Motor subtype changes in early Parkinson's disease," en, *Parkinsonism & Related Disorders*, vol. 43, pp. 67-72, Oct. 2017.</footnote>

<footnote n="12">¹²J. Jankovic, M. McDermott, J. Carter, et al., "Variable expression of Parkinson's disease: A base-line analysis of the DATATOP cohort. The Parkinson Study Group," eng, *Neurology*, vol. 40, no. 10, pp. 1529-1534, Oct. 1990.</footnote>

<footnote n="13">¹³R. B. Postuma, A. E. Lang, J. F. Gagnon, et al., "How does parkinsonism start? Prodromal parkinsonism motor changes in idiopathic REM sleep behaviour disorder," *Brain*, vol. 135, no. 6, pp. 1860-1870, Jun. 2012.</footnote>
progresses. In this case, some advocate for a continuum of variability to describe the differences between individuals. Hence the need to precisely identify the types of heterogeneity: a subtype-based heterogeneity, relying on an underlying causal factor for the differences, and the natural disease heterogeneity which leads to a continuous spectrum of observed variability.

**Early detection:** The last challenge with neurodegenerative diseases is to identify early markers of the disease. This is crucial in the enrollment of prodromal subjects into studies, in order to increase our understanding of the early development of the disease mechanisms. This also induces the possibility of early treatment, which is necessary in order to cure such diseases. Indeed, trying to revert the loss of neurons is an impossible task. The treatment has to act to prevent the neuronal loss, hence the administration time has to be at an early stage. But detecting early signs of a neurodegenerative disease is challenging, since many risk factors<fnref n="14" /> are not specific to one neurodegenerative disease: depression, constipation, or anxiety just to name a few ones. More specific markers can be detected with brain imaging or a cerebro-spinal fluid (CSF) sample, however such measurements are too costly and invasive to be performed routinely in order to detect early disease signs. More recently clinical research has striven to discover blood-based biomarkers (BBBM) of neurodegeneration<fnref n="15" /> which would be an adequate balance between sensitivity -which is lower than with imaging markers- and cost as well as invasiveness. Some more subtle signs are also studied, for instance the lost of the olfactive sense for AD<fnref n="16" />, or sleep troubles for PD<fnref n="17" />.

In order to address all of these challenges, high expectations were put upon data-driven models. Longitudinal cohorts have been built to follow the slow progression of those diseases, whose duration can span dozens of years. The databases are said to be longitudinal in the sense that individuals have repeated measurements at various points in time, which we will refer to as visits. This allows for the modelling of the dynamics of the disease, which is the goal of the field of disease progression modelling. The hope is that the accumulation of enough data from numerous centers with a wide array of measurements will help us better understand the disease, and ultimately lead to the discovery of the needed biomarkers for early detection and targets for therapeutic interventions.

Measurements in the cohorts include data from various nature: biomarkers such as BBBM, CSF or imaging biomarkers; cognitive assessments; patient questionnaires; motor tests; genetic sequencing... A special mention has to go to imaging, which has enormously benefited from the last decades developments of Magnetic Resonance Imaging (MRI) and Positron Emission Tomography (PET). Brain imaging has been at the root of the identification of abnormalities related to the neurodegenerative diseases. Famous large public research cohorts include the Alzheimer's Disease Neuroimaging Initiative (ADNI)<fnref n="18" />, and Parkinson's Progression Markers Initiative (PPMI)<fnref n="19" />.

<footnote n="14">¹⁴T. Nedelec, B. Couvy-Duchesne, F. Monnet, et al., "Identifying health conditions associated with Alzheimer's disease up to 15 years before diagnosis: An agnostic study of French and British health records," English, *The Lancet Digital Health*, vol. 4, no. 3, e169-e178, Mar. 2022.</footnote>

<footnote n="15">¹⁵C. E. Teunissen, I. M. W. Verberk, E. H. Thijssen, et al., "Blood-based biomarkers for Alzheimer's disease: Towards clinical implementation," en, *The Lancet Neurology*, vol. 21, no. 1, pp. 66-77, Jan. 2022.</footnote>

<footnote n="16">¹⁶C. Murphy, "Olfactory and other sensory impairments in Alzheimer disease," en, *Nature Reviews Neurology*, vol. 15, no. 1, pp. 11-24, Jan. 2019.</footnote>

<footnote n="17">¹⁷R. B. Postuma, A. E. Lang, J. F. Gagnon, et al., "How does parkinsonism start? Prodromal parkinsonism motor changes in idiopathic REM sleep behaviour disorder," *Brain*, vol. 135, no. 6, pp. 1860-1870, Jun. 2012.</footnote>

<footnote n="18">¹⁸R. C. Petersen, P. S. Aisen, L. A. Beckett, et al., "Alzheimer's Disease Neuroimaging Initiative (ADNI): Clinical characterization," en, *Neurology*, vol. 74, no. 3, pp. 201-209, Jan. 2010.</footnote>

<footnote n="19">¹⁹K. Marek, D. Jennings, S. Lasch, et al., "The Parkinson Progression Marker Initiative (PPMI)," en, *Progress in Neurobiology*, Biological Markers for Neurodegenerative Diseases, vol. 95, no. 4, pp. 629-635, Dec. 2011.</footnote>
The growth of these longitudinal cohorts has spurred the development of the disease progression modeling field. These digital models try to extract patterns of progression from the data. Their structure often allows to distinguish between two levels: global effects that describe the disease progression for the whole population, and individualized description of the disease to allow for subject-specific modelling. The goals of disease progression models are manifold.

**Diagnosis and detection:** One hope is that, by leveraging the vast amount of data collected from the longitudinal cohorts, the models are able to detect small signs of the disease onset ahead of the diagnosis. On the other hand, they can also be used to numerically determine thresholds for the diagnosis of the disease. Such models could be used to help the clinician with the diagnosis, and more realistically they could be implemented in the general practitioners offices to suggest when to redirect the patient to an expert neurologist consultation.

**Prognosis:** Disease progression models often provide individual predictions in time, thus allowing to forecast individuals in the future. This task is often used as the motivation for a disease progression model design, as predicting is more familiar to the Machine Learning community. This lead to the creation of prediction challenges for the progression of neurodegenerative diseases, such as the TADPOLE challenge<fnref n="20" />. Even though this competitive framework has sparked many interesting modelling attempts, performance should not be the sole aim of a model.

**Interpretability:** The model should provide insights into the disease typical evolution. Being able to exhibit parameters responsible for the individual specificity can be used to perform covariate analysis for instance. A common milestone is the design of what we call a disease timeline. It is well-known that age is not the proper reference axis to compare patients. Indeed some may start the disease earlier, others progress at a faster pace, or maybe both at the same time. Two patients at the same age might not be at the same disease stage. Hence the need to temporally align patients on a common disease timeline. Once this has been done, the variability left in the disease presentation can be referred to as spatial variability, as opposed to the time variability which accounts for realigning patients in time. A specific effort in our work has been put on the distinction between temporal and spatial effects.

Once a disease progression model has proven to reliably perform any of the tasks mentioned previously, the next step consists in using this model to achieve the ultimate goal of the research in neurodegenerative disease: find a cure. A first hope that is quickly dismissed is to find a potential target for a drug. This is highly unlikely, as the models are based on data that has been collected for the specific needs of the cohort. Data are thus recorded because of their already known, or hypothesized, link to the neurodegenerative process. A second potential use of disease progression models is to better identify when, how and to whom to administrate a drug. This is particularly useful in the era of personalized medicine. A final application context is the enhancement of clinical trials. Recently the US Food and Drugs Administration has approved the use of PROCOVA<fnref n="21" />, a

<footnote n="20">²⁰R. V. Marinescu, N. P. Oxtoby, A. L. Young, et al., "TADPOLE Challenge: Accurate Alzheimer's disease prediction through crowdsourced forecasting of future data," eng, *PRedictive Intelligence in MEdicine. PRIME (Workshop)*, vol. 11843, pp. 1-10, Oct. 2019.</footnote>

<footnote n="21">²¹D. Bertolini, K. Arnemann, D. Hall, et al., "Machine Learning Enables Smaller ALS Clinical Trials (P1-13.003)," en, *Neurology*, vol. 98, no. 18 Supplement, May 2022.</footnote>
method where the use of a prognostic disease progression score in an ANCOVA allows for a more powered trial, or equivalently allows for a smaller sample size.

<figure><img src="image_7.png" /><figcaption>Figure 1: Schematic representation of disease progression models. Data from various sources is recorded repeatedly in time for a patient. The aggregation of multiple individual visits in a longitudinal database allows to train disease progression models. These models can be used in multiple applications.</figcaption></figure>
Disease progression modelling

Unfolding from the hypothetical model from<fnref n="22" />, disease progression models have emerged to validate the hypothesis and provide a data-driven model as evidence. A specific class of models focuses on the ordering of the biomarkers abnormality progression. These models are called Event Based Models (EBM), and date back to<fnref n="23" />. Their application to neurodegenerative diseases uses cross-sectional data to determine the ordering of events<fnref n="24" />. When markers are continuous, events are defined by setting abnormality thresholds, allowing to include any type of data in the models. The EBM have been extended to account for the heterogeneity in the SUbtyping and STAge INference (SUSTAIN) algorithm<fnref n="25" />. However the sequence of events is only providing disease stages, and the disease timeline is not explicitly modelled.

Many models have focused on the time alignment as the cornerstone to estimate the disease progression. An early work proposed linear mapping to extract the disease timeline from individual observations<fnref n="26" />. This was later improved by the introduction of time-warps<fnref n="27" />. The construction of the common disease timeline by mapping patients' age onto a same axis through what is called a time reparametrization considers that individuals' observations are snapshots of a global common progression pattern, as in<fnref n="28" />. In this work the average trajectory was estimated without any prior shape hypothesis, only using the realigned individual snapshots. This common disease timeline can be seen as a continuous version of disease staging, a disease progression score<fnref n="29" />.

The main class of longitudinal models used in disease progression modelling has traditionally been the class of mixed-effect models<fnref n="30" />,<fnref n="31" />. This type of models has a hierarchical structure allowing to differentiate between population shared characteristics, or fixed effects, and individual variability, or random effects. Due to the non-linear nature of the progression of most biomarkers, the models belong to the non-linear mixed-effect models<fnref n="32" />. A family of template curves is chosen for the

<footnote n="22">²²C. R. Jack, D. S. Knopman, W. J. Jagust, et al., "Hypothetical model of dynamic biomarkers of the Alzheimer's pathological cascade," English, *The Lancet Neurology*, vol. 9, no. 1, pp. 119–128, Jan. 2010.</footnote>

<footnote n="23">²³L. A. Beckett, "Maximum Likelihood Estimation in Mallows's Model Using Partially Ranked Data," en, in *Probability Models and Statistical Analyses for Ranking Data*, M. A. Fligner and J. S. Verducci, Eds., ser. Lecture Notes in Statistics, New York, NY: Springer, 1993, pp. 92–107.</footnote>

<footnote n="24">²⁴H. M. Fonteijn, M. Modat, M. J. Clarkson, et al., "An event-based model for disease progression and its application in familial Alzheimer's disease and Huntington's disease," eng, *NeuroImage*, vol. 60, no. 3, pp. 1880–1889, Apr. 2012.</footnote>

<footnote n="25">²⁵A. L. Young, R. V. Marinescu, N. P. Oxtoby, et al., "Uncovering the heterogeneity and temporal complexity of neurodegenerative diseases with Subtype and Stage Inference," en, *Nature Communications*, vol. 9, no. 1, pp. 1–16, Oct. 2018.</footnote>

<footnote n="26">²⁶E. Yang, M. Farnum, V. Lobanov, et al., "Quantifying the pathophysiological timeline of Alzheimer's disease," eng, *Journal of Alzheimer's disease: JAD*, vol. 26, no. 4, pp. 745–753, 2011.</footnote>

<footnote n="27">²⁷S. Durrleman, X. Pennec, A. Trouvé, et al., "Toward a Comprehensive Framework for the Spatiotemporal Statistical Analysis of Longitudinal Shape Data," en, *International Journal of Computer Vision*, vol. 103, no. 1, pp. 22–59, May 2013.</footnote>

<footnote n="28">²⁸M. C. Donohue, H. Jacqmin-Gadda, M. Le Goff, et al., "Estimating long-term multivariate progression from short-term data," *Alzheimer's & dementia: the journal of the Alzheimer's Association*, vol. 10, no. 0, S400–S410, Oct. 2014.</footnote>

<footnote n="29">²⁹V. Kmetzsch, E. Becker, D. Saracino, et al., "Disease Progression Score Estimation From Multimodal Imaging and MicroRNA Data Using Supervised Variational Autoencoders," *IEEE Journal of Biomedical and Health Informatics*, vol. 26, no. 12, pp. 6024–6035, Dec. 2022.</footnote>

<footnote n="30">³⁰N. M. Laird and J. H. Ware, "Random-effects models for longitudinal data," eng, *Biometrics*, vol. 38, no. 4, pp. 963–974, Dec. 1982.</footnote>

<footnote n="31">³¹M. J. Lindstrom and D. M. Bates, "Newton—Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data," *Journal of the American Statistical Association*, vol. 83, no. 404, pp. 1014–1022, Dec. 1988.</footnote>

<footnote n="32">³²J. Serroyen, G. Molenberghs, G. Verbeke, et al., "Nonlinear Models for Longitudinal Data," *The American*</footnote>
progression, for instance logistic curves, and with an affine time-reparametrization the individuals
are mapped onto the average trajectory, as in<fnref n="33" />. Other choices of template curves are possible, as
the ones proposed by<fnref n="34" />. In the R-package for latent classes mixed effect models (LCMM)<fnref n="35" />, the
non-linear model can be combined with the joint estimation of subtypes. Non-linear mixed effect
models describing the progression of continuous markers can be coupled with event modelling, in
what is called joint models<fnref n="36" />. An interesting approach used Gaussian Processes (GP) to account for
the random effects of the mixed effect model<fnref n="37" />,<fnref n="38" />. The Disease Course Mapping (DCM) approach
taken in<fnref n="39" /> relies on a geometric model assuming that the observations lie on a Riemannian manifold.
It is inspired from a prior geometric model<fnref n="40" />, which was computationally intractable because the
noise was defined in the manifold.

These non-linear models have also been extended to high-dimensional modalities such as imag-
ing. Extending<fnref n="41" />, the authors in<fnref n="42" /> applied the model to images in ADNI. DIVE<fnref n="43" /> is voxel-based
approach for modelling the evolution of the brain in neurodegenerative diseases. The GP model
was extended to shape modelling in<fnref n="44" />. In<fnref n="45" />, the DCM model was applied to brain imaging, where
the brain is modelled as a mesh and embedded in a shape manifold.

Modelling the dynamic progression of the biomarkers has also been performed using dynamical
models. They are based on Ordinary Differential Equations (ODE)<fnref n="46" />,<fnref n="47" /> or discretized versions of

<footnote>*Statistician*, vol. 63, no. 4, pp. 378-388, Nov. 2009.</footnote>

<footnote n="33">³³B. M. Jedynak, A. Lang, B. Liu, *et al.*, "A Computational Neurodegenerative Disease Progression Score: Method and Results with the Alzheimer's Disease Neuroimaging Initiative Cohort," *NeuroImage*, vol. 63, no. 3, pp. 1478-1486, Nov. 2012.</footnote>

<footnote n="34">³⁴L. L. Raket, “Statistical Disease Progression Modeling in Alzheimer Disease,” *Frontiers in Big Data*, vol. 3, 2020.</footnote>

<footnote n="35">³⁵C. Proust-Lima, V. Philipps, and B. Liquet, “Estimation of Extended Mixed Models Using Latent Classes and Latent Processes: The R Package lcm,” en, *Journal of Statistical Software*, vol. 78, pp. 1–56, Jun. 2017.</footnote>

<footnote n="36">³⁶B. He and S. Luo, “Joint modeling of multivariate longitudinal measurements and survival data with applications to Parkinson’s disease,” eng, *Statistical Methods in Medical Research*, vol. 25, no. 4, pp. 1346–1358, Aug. 2016.</footnote>

<footnote n="37">³⁷M. Lorenzi, X. Pennec, G. B. Frisoni, et al., “Disentangling normal aging from Alzheimer’s disease in structural magnetic resonance images,” en, *Neurobiology of Aging*, Novel Imaging Biomarkers for Alzheimer’s Disease and Related Disorders (NIBAD), vol. 36, S42–S52, Jan. 2015.</footnote>

<footnote n="38">³⁸M. Lorenzi, M. Filippone, G. B. Frisoni, et al., “Probabilistic disease progression modeling to characterize diagnostic uncertainty: Application to staging and prediction in Alzheimer’s disease,” en, *NeuroImage*, Mapping diseased brains, vol. 190, pp. 56–68, Apr. 2019.</footnote>

<footnote n="39">³⁹J.-B. Schiratti, S. Allassonnière, A. Routier, et al., “A Mixed-Effects Model with Time Reparametrization for Longitudinal Univariate Manifold-Valued Data,” en, in *Information Processing in Medical Imaging*, S. Ourselin, D. C. Alexander, C.-F. Westin, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2015, pp. 564–575.</footnote>

<footnote n="40">⁴⁰T. Fletcher, “Geodesic Regression on Riemannian Manifolds,” en, Sep. 2011, p. 75.</footnote>

<footnote n="41">⁴¹B. M. Jedynak, A. Lang, B. Liu, et al., "A Computational Neurodegenerative Disease Progression Score: Method and Results with the Alzheimer's Disease Neuroimaging Initiative Cohort," *NeuroImage*, vol. 63, no. 3, pp. 1478-1486, Nov. 2012.</footnote>

<footnote n="42">⁴²M. Bilgel, J. L. Prince, D. F. Wong, et al., “A multivariate nonlinear mixed effects model for longitudinal image analysis: Application to amyloid imaging,” en, *NeuroImage*, vol. 134, pp. 658–670, Jul. 2016.</footnote>

<footnote n="43">⁴³R. V. Marinescu, A. Eshaghi, M. Lorenzi, et al., “DIVE: A spatiotemporal progression model of brain pathology in neurodegenerative disorders,” en, *NeuroImage*, vol. 192, pp. 166–177, May 2019.</footnote>

<footnote n="44">⁴⁴C. Abi Nader, N. Ayache, P. Robert, et al., “Monotonic Gaussian Process for spatio-temporal disease progression modeling in brain imaging data,” eng, *NeuroImage*, vol. 205, p. 116266, Jan. 2020.</footnote>

<footnote n="45">⁴⁵I. Koval, “Learning Multimodal Digital Models of Disease Progression from Longitudinal Data: Methods & Algorithms for the Description, Prediction and Simulation of Alzheimer’s Disease Progression,” en, Ph.D. dissertation, Institut Polytechnique de Paris, Jan. 2020.</footnote>

<footnote n="46">⁴⁶N. P. Oxtoby, A. L. Young, D. M. Cash, et al., “Data-driven models of dominantly-inherited Alzheimer’s disease progression,” *Brain*, vol. 141, no. 5, pp. 1529–1544, May 2018.</footnote>

<footnote n="47">⁴⁷K. Lahouel, M. Wells, V. Rielly, et al., *Learning nonparametric ordinary differential equations from noisy data*,</footnote>
ODE<sup>48</sup>, where the field of the ODE is estimated from the whole cohort and is akin to the fixed effects of the mixed models while the initial conditions of an individual patients trajectory account for the random effects. Neural ODE<sup>49</sup> use a Deep Learning network to learn the field. Other Deep Learning approaches using Recurrent Neural Networks have also been applied for disease progression<sup>50</sup>. These Deep Learning methods are often less interpretable. This is quite the opposite to the multifactorial causal model proposed in<sup>51</sup>. This causal model allows to perform interventions, which can be leveraged to emulate therapeutic intervention.

All those statistical methods have in common the ability to describe the progression of the disease at a population level. Most of them can also compute an individualized trajectory for every patient. The main differences lie in their interpretability. Models with an explicit time-reparametrization provide a separation of temporal and spatial variability, which is an asset to us. Semi-parametric or non-parametric models, such as the dynamical models based on ODE or Deep Learning based methods, often lack on the interpretability side. Our work focuses on the geometrical approach of the DCM, which provides both the spatio-temporal disentanglement and meaningful parameters.

Parkinson's disease

Most of the work in this thesis has been developed for applications to Parkinson's disease data, so we will provide a brief overview of this disease, its hypothesized causes, its manifestation and its treatments.

PD was first described by James Parkinson in 1817. This neurodegenerative disease is the second most prevalent after AD. Its prevalence increases drastically with age, from being almost null before 50 years to close to 2% after 80 years old<sup>52</sup>. Men are twice more likely to develop the disease compared to women.

PD is characterized by the loss of dopaminergic neurons, starting in the substantia nigra, a region in the brainstem containing almost exclusively dopaminergic neurons. Dopamine is a neurotransmitter involved in the reward system, but is also responsible for the activation of the motor system. The decrease in dopamine levels due to neurodegeneration entails all the motor symptoms of PD that we are familiar with: bradykinesia (slowness of movement), muscle rigidity, tremors, which are the main symptoms<sup>53</sup>, but many others are also observed including falls, slowness of speech... The loss of dopaminergic neurons spreads in time to other regions of the brain, inducing non-motor symptoms, such as apathy, hallucinations or depression. However the onset of symptom is a very late landmark of neuronal death. When PD is diagnosed, already 80% of the dopaminergic neurons in the substantia nigra are already lost. The neurodegeneration is hypothesized to come

<footnote>Feb. 2023.</footnote>

<footnote n="48"><sup>48</sup>B. O. Taddé, H. Jacqmin-Gadda, J.-F. Dartigues, et al., “Dynamic modeling of multivariate dimensions and their temporal relationships using latent processes: Application to Alzheimer’s disease,” en, *Biometrics*, vol. 76, no. 3, pp. 886–899, 2020.</footnote>

<footnote n="49"><sup>49</sup>R. T. Q. Chen, Y. Rubanova, J. Bettencourt, et al., *Neural Ordinary Differential Equations*, Dec. 2019.</footnote>

<footnote n="50"><sup>50</sup>M. Nguyen, T. He, L. An, et al., “Predicting Alzheimer’s disease progression using deep recurrent neural networks,” en, *NeuroImage*, vol. 222, p. 117–203, Nov. 2020.</footnote>

<footnote n="51"><sup>51</sup>Y. Iturria-Medina, F. M. Carbonell, R. C. Sotero, et al., “Multifactorial causal model of brain (dis)organization and therapeutic intervention: Application to Alzheimer’s disease,” en, *NeuroImage*, vol. 152, pp. 60–77, May 2017.</footnote>

<footnote n="52"><sup>52</sup>A. Elbaz, L. Carcaillon, S. Kab, et al., “Epidemiology of Parkinson’s disease,” en, *Revue Neurologique*, Neuropidemiology, vol. 172, no. 1, pp. 14–26, Jan. 2016.</footnote>

<footnote n="53"><sup>53</sup>J. Jankovic, M. McDermott, J. Carter, et al., “Variable expression of Parkinson’s disease: A base-line analysis of the DATATOP cohort. The Parkinson Study Group,” eng, *Neurology*, vol. 40, no. 10, pp. 1529–1534, Oct. 1990.</footnote>
from the accumulation of the α-synuclein protein, which is currently the main target for drug de-
velopment<sup>54</sup>. The loss of dopaminergic neurons can be observed with dopamine transporter Single
Photon Emission Computed Tomography (SPECT) scan.

As mentioned in the challenges about neurodegenerative diseases, the definition of PD is com-
plex. It is currently clinical, at the neurologist discretion, but there is no clear-cut guideline for
diagnosis. The symptoms are similar to Lewy-body disease, which leads to the early diagnosis
being revised sometimes. However we are still far from a biological definition of the disease. The
heterogeneity of the disease plays an important role, starting from the symptom onset<sup>55</sup>. A first
distinction can be made based on the location of the first motor signs. 60% of the patients are
first impaired on one side of the body, the symptoms spreading slowly towards the rest of the
body in time. Attempts at subtyping have tried to distinguish early prodromal subtypes, based
on sleep behaviours<sup>56</sup>. One of the most common classification based on motor symptoms is the
Tremor Dominant (TD) vs Postural Instability and Gait Disorder (PIGD) system proposed in<sup>57</sup>.
However these subtypes do not encompass all of the patients, and they are not stable in time<sup>58,59</sup>.
Some particular subtypes can be derived from genetic variants GBA and LRKK2, which account
for about 5% of the total PD cases, the others being idiopathic.

PD has benefited from the discovery of a treatment in the early 50's, called Levodopa. This is a dopamine precursor, and when taken regularly it makes up for the lack of dopamine production in the brain. Other drugs based on dopaminergic input have since been developed such as dopamin agonists, or inhibitors of recapture including IMAOB and ICOMT. They all rely on the same mechanism, with different side-effects. Levodopa is for instance known to cause dyskinesia, i.e. spontaneous movements, and has also been related to freezing in some studies<sup>60</sup>. Dopamin agonists are more likely to cause Impulse Control Disorders (ICD), which causes addictive behaviours such as gambling. The treatment for PD patients nowadays is a cocktail of several dopaminergic drugs, experts are still debating what is the best treatment policy<sup>61</sup>. Many factors can influence the answer of a patient to the treatment, and its likelihood of developing negative side-effects.

The dopaminergic treatment is purely symptomatic however, it has not shown any improvement
in the long run<sup>62</sup>. Its efficacy diminishes as the disease progresses towards later stages, leading to
fluctuations between "ON" and "OFF" states. "ON" states correspond to the normally functioning
drug state, where motor symptoms are alleviated. "OFF" states on the other hand occur when the
treatment effect wears out. In later stages of PD the half-life of the dopaminergic drug becomes

<footnote n="54">⁵⁴C. R. Fields, N. Bengoa-Vergniory, and R. Wade-Martins, "Targeting Alpha-Synuclein as a Therapy for Parkinson's Disease," *Frontiers in Molecular Neuroscience*, vol. 12, 2019.</footnote>

<footnote n="55">⁵⁵W. J. Zetusky, J. Jankovic, and F. J. Pirozzolo, "The heterogeneity of Parkinson's disease: Clinical and prognostic implications," eng, *Neurology*, vol. 35, no. 4, pp. 522–526, Apr. 1985.</footnote>

<footnote n="56">⁵⁶R. B. Postuma, A. E. Lang, J. F. Gagnon, et al., "How does parkinsonism start? Prodromal parkinsonism motor changes in idiopathic REM sleep behaviour disorder," *Brain*, vol. 135, no. 6, pp. 1860–1870, Jun. 2012.</footnote>

<footnote n="57">⁵⁷J. Jankovic, M. McDermott, J. Carter, et al., "Variable expression of Parkinson's disease: A base-line analysis of the DATATOP cohort. The Parkinson Study Group," eng, *Neurology*, vol. 40, no. 10, pp. 1529–1534, Oct. 1990.</footnote>

<footnote n="58">⁵⁸T. Simuni, C. Caspell-Garcia, C. Coffey, et al., "How stable are Parkinson's disease subtypes in de novo patients: Analysis of the PPMI cohort?" en, *Parkinsonism & Related Disorders*, vol. 28, pp. 62–67, Jul. 2016.</footnote>

<footnote n="59">⁵⁹R. S. Eisinger, C. W. Hess, D. Martinez-Ramirez, et al., "Motor subtype changes in early Parkinson's disease," en, *Parkinsonism & Related Disorders*, vol. 43, pp. 67–72, Oct. 2017.</footnote>

<footnote n="60">⁶⁰M. Gilat, N. D'Cruz, P. Ginis, et al., "Freezing of gait and levodopa," English, *The Lancet Neurology*, vol. 20, no. 7, pp. 505–506, Jul. 2021.</footnote>

<footnote n="61">⁶¹R. M. A. de Bie, C. E. Clarke, A. J. Espay, et al., "Initiation of pharmacological therapy in Parkinson's disease: When, why, and how," eng, *The Lancet. Neurology*, vol. 19, no. 5, pp. 452–461, May 2020.</footnote>

<footnote n="62">⁶²R. Cilia, E. Cereda, A. Akpalu, et al., "Natural history of motor symptoms in Parkinson's disease and the long-duration response to levodopa," *Brain*, vol. 143, no. 8, pp. 2490–2501, Aug. 2020.</footnote>
really short due to the quasi-inexisting production by the brain.

Another treatment called Deep Brain Stimulation implants an electrode in the amygdala to directly stimulate the brain. This treatment has shown to be a potent alternative to dopaminergic therapy<sup>63</sup>, but still remains quite rare due to its invasiveness. The dopaminergic treatment has long overshadowed the research for a disease-modifying drug, therefore the research field pales in comparison with its counterpart in AD.

Longitudinal cohorts used in this thesis work include PPMI<sup>64</sup> and a french cohort from the NS-PARK consortium. The data in PD cohorts is mostly composed of cognitive and motor assessments. Historically the Hoehn and Yahr scale<sup>65</sup> (from 0 to 5) was used to describe the disease in successive stages. Today the most frequently used scale for disease progression is the Movement Disorder Society Unified Parkinson's Disease Rating Scale (MDS-UPDRS)<sup>66</sup>, which is an aggregation of many items rating multiple aspects of the disease, including a motor examination and an assessment of the disease impact on the daily life of the patient.

## Thesis goals

The work in this thesis was built upon the Disease Course Mapping model introduced in<sup>67</sup>. We have described some key challenges raised by the neurodegenerative diseases and we will extend the DCM model to address those challenges. The DCM model has already been built to account for time alignment and spatio-temporal decoupling. However its design is far from perfect. First it is limited to continuous data. Moreover, the DCM model does not account for the heterogeneity of disease subtypes. Finally it is not geared for applications in drug development. Our brief overview of Parkinson's disease contrasts with those limitations. Symptoms, staging and rating scales are ubiquitous, and these are non-continuous variables. Subtypes are paramount in modelling the variability of PD. Finally all the data in PD cohorts comes from patients under dopaminergic treatment, thus requiring a specific model for it. We will therefore provide the following extensions of the DCM model:

* applicability to discrete data, including binary and ordinal observations

* mixture structure to uncover subtypes and better explain disease heterogeneity

* treatment models linked to the DCM to account for treatment effect and pave the way for the development of disease-modifying drugs

<footnote n="63"><sup>63</sup>A. L. Benabid, "Deep brain stimulation for Parkinson's disease," en, *Current Opinion in Neurobiology*, vol. 13, no. 6, pp. 696-706, Dec. 2003.</footnote>

<footnote n="64"><sup>64</sup>K. Marek, D. Jennings, S. Lasch, et al., "The Parkinson Progression Marker Initiative (PPMI)," en, *Progress in Neurobiology*, Biological Markers for Neurodegenerative Diseases, vol. 95, no. 4, pp. 629-635, Dec. 2011.</footnote>

<footnote n="65"><sup>65</sup>M. M. Hoehn and M. D. Yahr, "Parkinsonism: Onset, progression, and mortality," en, *Neurology*, vol. 17, no. 5, pp. 427-427, May 1967.</footnote>

<footnote n="66"><sup>66</sup>C. G. Goetz, B. C. Tilley, S. R. Shaftman, et al., "Movement Disorder Society-sponsored revision of the Unified Parkinson's Disease Rating Scale (MDS-UPDRS): Scale presentation and clinimetric testing results: MDS-UPDRS: Clinimetric Assessment," en, *Movement Disorders*, vol. 23, no. 15, pp. 2129-2170, Nov. 2008.</footnote>

<footnote n="67"><sup>67</sup>J.-B. Schiratti, S. Allassonnière, A. Routier, et al., "A Mixed-Effects Model with Time Reparametrization for Longitudinal Univariate Manifold-Valued Data," en, in *Information Processing in Medical Imaging*, S. Ourselin, D. C. Alexander, C.-F. Westin, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2015, pp. 564-575.</footnote>
All our methodological contributions are motivated by practical needs encountered in real-life data. They are thus accompanied by applications to neurodegenerative cohorts. We used the interpretability of the model to provide clinical insights in the analysis of our experimental results.

Manuscript overview

The first chapter 1 introduces the Disease Course Mapping model and discusses its limitations. The second part of the chapter presents Geodesic Bending, an extension of the DCM model which introduces more flexibility in the model by learning the shape of the typical trajectory of the disease.

The second chapter 2 presents our extension of the DCM model to binary and ordinal data. We provide in-depth applications to Parkinson's disease cohorts, as well as a parallel with the Item Response Theory.

The third chapter 3 addresses the heterogeneity challenge of neurodegenerative diseases. Our contribution takes the form of two mixture variants for the DCM model. We showcase the differences between the two with an application to Alzheimer's disease.

The fourth chapter 4 focuses on treatment modelling. We introduce two modelling approaches. The first one consists in a generic framework using any disease progression model to estimate treatment effect from patients only, without requiring data from controls nor placebo. The second approach is called Piece-wise Geodesic model and is targeted at disease-modifying drugs, providing a framework to learn changes in disease trajectory due to the treatment.

The fifth chapter 5 summarizes all our contributions to the estimation algorithm across all of the previous extensions. It presents the MCMC-SAEM algorithm used for the DCM estimation. We then introduce our contributions in the form of specific variants of this algorithm.
# Introduction

## Motivation

À la suite des progrès médicaux et sanitaires réalisés au cours des derniers siècles, nos sociétés ont considérablement gagné en espérance de vie. Avec l'émergence d'une population de plus en plus âgée, de nouveaux défis en matière de santé ont surgi, notamment le cancer et les maladies neurodégénératives. En se concentrant sur ces dernières, cette catégorie de pathologies englobe un large spectre de maladies, se caractérisant par un dysfonctionnement du système nerveux périphérique ou central. Actuellement, les maladies neurodégénératives sont la deuxième cause de décès dans le monde<fnref n="68" />. On estime qu'une femme sur deux et un homme sur trois seront touchés au cours de leur vie. La démence affecte une personne sur 14 après l'âge de 65 ans et une personne sur 6 après 80 ans, la maladie d'Alzheimer étant la cause la plus fréquente de démence, représentant environ 70% de tous les cas de démence. La deuxième maladie neurodégénérative la plus courante est la maladie de Parkinson, qui compte plus de 12 millions de patients en Europe. Ces chiffres sont appelés à augmenter dans les années à venir, car les prévisions démographiques prévoient une augmentation de la population de 65 ans et plus de 18% à 24% en 2060 (statistiques de la population française<fnref n="69" />). Cependant, ces maladies neurodégénératives soulèvent plusieurs problèmes dans leur étude.

**Comprendre les maladies:** Contrairement aux maladies infectieuses, où la cause est facilement identifiable et les mécanismes pathogènes peuvent être étudiés, les maladies neurodégénératives restent largement méconnues. Il n'y a ni bactérie ni virus à blâmer, aucun début clair tel qu'une infection, aucune cure ni vaccin. Leur physiopathologie repose sur des symptômes et des évaluations cliniques des capacités cognitives ou motrices, qui résultent d'une neurodégénérescence tardive. Cependant, le processus de neurodégénérescence débute bien plus tôt, passant inaperçu jusqu'à ce qu'il soit trop tard. Par exemple, dans la maladie d'Alzheimer ou la maladie de Parkinson, la maladie peut durer des décennies, et des études montrent que les mécanismes biologiques responsables des dommages neuronaux ont commencé plus de vingt ans avant le diagnostic<fnref n="70" />. Cette évolution cachée de la neurodégénérescence est la principale raison pour laquelle la recherche médicale a du retard dans la compréhension de ces maladies. Les patients diagnostiqués ne représentent qu'une phase tardive de la maladie, lorsque la plupart des dommages ont déjà été causés. L'étude de l'histoire naturelle de la maladie est donc entravée par un manque de sujets présentant un stade

<footnote n="68">⁶⁸J. Dumurgier and C. Tzourio, "Epidemiology of neurological diseases in older adults," en, *Revue Neurologique*, vol. 176, no. 9, pp. 642-648, Nov. 2020.</footnote>

<footnote n="69">⁶⁹*Projections de population à l'horizon 2060 - Insee Première - 1320.*</footnote>

<footnote n="70">⁷⁰R. B. Postuma, A. E. Lang, J. F. Gagnon, et al., "How does parkinsonism start? Prodromal parkinsonism motor changes in idiopathic REM sleep behaviour disorder," *Brain*, vol. 135, no. 6, pp. 1860-1870, Jun. 2012.</footnote>
précoce de la maladie, dans un état appelé prodromique. Le recrutement de personnes prodromiques est particulièrement difficile en raison d'une autre lacune dans nos connaissances sur les maladies neurodégénératives : les causes restent en grande partie inconnues. Des facteurs de risque ont été identifiés<sup>71,72</sup>, y compris des facteurs environnementaux ou génétiques. Cependant, ils ne représentent qu'une partie des cas de la maladie, et leur présence ne signifie pas automatiquement que l'individu à risque développera la pathologie, d'où la nuance entre un facteur de risque et une cause. Nous sommes souvent confrontés à la grande majorité de patients qui sont affectés spontanément, ce que l'on appelle la forme idiopathique de la maladie.

L'étude de ces maladies a conduit à des modèles hypothétiques tentant de façonner notre compréhension de leur progression. Dans la maladie d'Alzheimer, le célèbre mécanisme de cascade amyloïde a été proposé<sup>73</sup>. Ce modèle mettait en évidence le rôle de la protéine β-amyloïde, dont l'accumulation progressive dans le cerveau entraînait la mort neuronale, et finalement la démence. Ce modèle hypothétique a façonné la recherche dans le domaine depuis lors, au point où les protéines amyloïdes agrégées sont la cible principale des essais thérapeutiques. L'étude de cette cascade sur les aspects observables de la maladie a également conduit à un autre modèle hypothétique bien connu dans la maladie d'Alzheimer. Introduit dans une étude<sup>74</sup>, les auteurs ont proposé que l'évolution des marqueurs de la maladie, tels que la charge amyloïde, les marqueurs de la mort neuronale et le déclin cognitif, suivrait une progression logistique du normal à l'anormal. Des modèles similaires ont également été proposés dans la maladie de Parkinson, notamment les stades de Braak<sup>75</sup> et plus récemment l'hypothèse du cerveau d'abord et du corps d'abord<sup>76</sup>. La partie cruciale de ces modèles descriptifs consiste à fournir un *ordre* de dégradation des marqueurs, qui est ensuite utilisé pour définir les stades de la maladie. La stratification en stades est essentielle, comme nous le verrons lors de la discussion sur le développement d'un traitement.

**Définir la maladie:** Cette compréhension naissante des processus biologiques sous-jacents à la physiopathologie observable de la maladie remet en question même la définition de ces maladies. Jusqu'à récemment, la notion de maladie était liée aux manifestations cliniques, c'est-à-dire aux symptômes, et le diagnostic était évalué par un neurologue. Cependant, le patient avait probablement la maladie avant l'apparition des symptômes. Par conséquent, les développements récents plaident en faveur d'une définition biologique des maladies, basée sur des mesures objectives de biomarqueurs qui permettraient un diagnostic préalable au diagnostic clinique. Dans la maladie d'Alzheimer, une telle définition repose sur la présence d'amyloïde dans le cerveau, ou de la charge en tau - une autre protéine dont le mauvais repliement induit la formation de noeuds et de fibrilles toxiques pour les neurones - et les patients sont ainsi décrits selon une classification A/T/N (pour

<footnote n="71"><sup>71</sup>G. Livingston, J. Huntley, A. Sommerlad, et al., "Dementia prevention, intervention, and care: 2020 report of the Lancet Commission," English, *The Lancet*, vol. 396, no. 10248, pp. 413-446, Aug. 2020.</footnote>

<footnote n="72"><sup>72</sup>T. Nedelec, B. Couvy-Duchesne, F. Monnet, et al., "Identifying health conditions associated with Alzheimer's disease up to 15 years before diagnosis: An agnostic study of French and British health records," English, *The Lancet Digital Health*, vol. 4, no. 3, e169-e178, Mar. 2022.</footnote>

<footnote n="73"><sup>73</sup>J. A. Hardy and G. A. Higgins, "Alzheimer's disease: The amyloid cascade hypothesis," eng, *Science (New York, N.Y.)*, vol. 256, no. 5054, pp. 184-185, Apr. 1992.</footnote>

<footnote n="74"><sup>74</sup>C. R. Jack, D. S. Knopman, W. J. Jagust, et al., "Hypothetical model of dynamic biomarkers of the Alzheimer's pathological cascade," English, *The Lancet Neurology*, vol. 9, no. 1, pp. 119-128, Jan. 2010.</footnote>

<footnote n="75"><sup>75</sup>H. Braak, K. D. Tredici, U. Rüb, et al., "Staging of brain pathology related to sporadic Parkinson's disease," en, *Neurobiology of Aging*, vol. 24, no. 2, pp. 197-211, Mar. 2003.</footnote>

<footnote n="76"><sup>76</sup>J. Horsager, K. B. Andersen, K. Knudsen, et al., "Brain-first versus body-first Parkinson's disease: A multi-modal imaging case-control study," *Brain*, vol. 143, no. 10, pp. 3077-3088, Oct. 2020.</footnote>
Amyloïde/Tau/Neurodégénération)$^{77}$.

**Décrire l'hétérogénéité:** Cet dernier exemple de classification chez les patients atteints de la maladie d'Alzheimer met en évidence la grande variété de manifestations possibles de la maladie. Souvent, une maladie neurodégénérative englobe divers sous-types, parfois en raison de variants génétiques spécifiques. Cette hétérogénéité dans le spectre des maladies explique également pourquoi il est difficile de définir correctement la maladie. Dans certains cas, l'hétérogénéité est inévitable car l'étude de sous-groupes plus restreints de patients laisserait trop peu de données pour étudier la maladie et trop peu de sujets potentiels pour être recrutés dans un essai clinique. Il est donc toujours nécessaire de tenir compte de la part commune de la physiologie de la maladie ainsi que des différences entre les sous-types. La distinction entre les sous-types n'est pas toujours claire, car il n'y a pas toujours de variant génétique ou de marqueur distinctif. De nombreuses tentatives ont été faites pour la sous-typage de la maladie de Parkinson, par exemple$^{78,79,80}$, mais les frontières entre les sous-groupes sont floues et pas toujours stables au fur et à mesure de la progression de la maladie. Dans ce cas, certains plaident en faveur d'un continuum de variabilité pour décrire les différences entre les individus. D'où la nécessité d'identifier précisément les types d'hétérogénéité : une hétérogénéité basée sur les sous-types, reposant sur un facteur causal sous-jacent pour les différences, et l'hétérogénéité naturelle de la maladie qui conduit à un spectre continu de variabilité observée.

**Détection précoce:** Le dernier défi des maladies neurodégénératives consiste à identifier des marqueurs précoces de la maladie. Cela est crucial pour l'inclusion de sujets prodromiques dans les études, afin d'accroître notre compréhension du développement précoce des mécanismes de la maladie. Cela ouvre également la possibilité d'un traitement précoce, qui est nécessaire pour guérir de telles maladies. En effet, tenter de rétablir la perte de neurones est une tâche impossible. Le traitement doit agir pour prévenir la perte neuronale, d'où la nécessité d'une administration précoce. Cependant, détecter les premiers signes d'une maladie neurodégénérative est un défi, car de nombreux facteurs de risque$^{81}$ ne sont pas spécifiques à une maladie neurodégénérative particulière : la dépression, la constipation ou l'anxiété, pour n'en nommer que quelques-uns. Des marqueurs plus spécifiques peuvent être détectés par imagerie cérébrale ou par prélèvement de liquide céphalorachidien (LCR), cependant de telles mesures sont trop coûteuses et invasives pour être réalisées systématiquement afin de détecter les premiers signes de la maladie. Plus récemment, la recherche clinique s'est efforcée de découvrir des biomarqueurs de neurodégénérescence basés sur le sang (BBBM)$^{82}$, ce qui représenterait un équilibre adéquat entre sensibilité - qui est inférieure à celle des marqueurs d'imagerie - et coût ainsi que l'aspect invasif. Certains signes plus subtils sont

<footnote>$^{77}$C. R. Jack, D. A. Bennett, K. Blennow, et al., "A/T/N: An unbiased descriptive classification scheme for Alzheimer disease biomarkers," eng, *Neurology*, vol. 87, no. 5, pp. 539-547, Aug. 2016.</footnote>

<footnote>$^{78}$R. S. Eisinger, C. W. Hess, D. Martinez-Ramirez, et al., "Motor subtype changes in early Parkinson's disease," en, *Parkinsonism & Related Disorders*, vol. 43, pp. 67-72, Oct. 2017.</footnote>

<footnote>$^{79}$J. Jankovic, M. McDermott, J. Carter, et al., "Variable expression of Parkinson's disease: A base-line analysis of the DATATOP cohort. The Parkinson Study Group," eng, *Neurology*, vol. 40, no. 10, pp. 1529-1534, Oct. 1990.</footnote>

<footnote>$^{80}$R. B. Postuma, A. E. Lang, J. F. Gagnon, et al., "How does parkinsonism start? Prodromal parkinsonism motor changes in idiopathic REM sleep behaviour disorder," *Brain*, vol. 135, no. 6, pp. 1860-1870, Jun. 2012.</footnote>

<footnote>$^{81}$T. Nedelec, B. Couvy-Duchesne, F. Monnet, et al., "Identifying health conditions associated with Alzheimer's disease up to 15 years before diagnosis: An agnostic study of French and British health records," English, *The Lancet Digital Health*, vol. 4, no. 3, e169-e178, Mar. 2022.</footnote>

<footnote>$^{82}$C. E. Teunissen, I. M. W. Verberk, E. H. Thijssen, et al., "Blood-based biomarkers for Alzheimer's disease: Towards clinical implementation," en, *The Lancet Neurology*, vol. 21, no. 1, pp. 66-77, Jan. 2022.</footnote>
également étudiés, par exemple la perte de l'odorat dans la maladie d'Alzheimer<sup>83</sup> ou les troubles du sommeil dans la maladie de Parkinson<sup>84</sup>.

Pour relever tous ces défis, de grandes attentes ont été placées dans les modèles basés sur les données. Des cohortes longitudinales ont été constituées pour suivre la lente progression de ces maladies, dont la durée peut s'étendre sur des dizaines d'années. Les bases de données sont qualifiées de longitudinales dans le sens où les individus ont des mesures répétées à différents moments, que nous appellerons visites. Cela permet de modéliser la dynamique de la maladie, qui est l'objectif du domaine de la modélisation de la progression des maladies. L'espoir est que l'accumulation de suffisamment de données provenant de nombreux centres et comportant une large gamme de mesures nous aidera à mieux comprendre la maladie et à découvrir en fin de compte les biomarqueurs nécessaires pour une détection précoce et les cibles pour des interventions thérapeutiques.

Les mesures dans les cohortes comprennent des données de diverses natures : des biomarqueurs tels que BBBM, LCR ou des biomarqueurs d'imagerie ; des évaluations cognitives ; des questionnaires pour les patients ; des tests moteurs ; des séquençages génétiques... Une mention spéciale doit être faite pour l'imagerie, qui a énormément bénéficié des développements des dernières décennies en imagerie par résonance magnétique (IRM) et en tomographie par émission de positrons (TEP). L'imagerie cérébrale a été à l'origine de l'identification des anomalies liées aux maladies neurologique de la maladie d'Alzheimer (ADNI)<sup>85</sup> et l'initiative pour les marqueurs de progression de la maladie de Parkinson (PPMI)<sup>86</sup>.

L'essor de ces cohortes longitudinales a stimulé le développement du domaine de la modélisation de la progression des maladies. Ces modèles numériques tentent d'extraire des schémas de progression à partir des données. Leur structure permet souvent de distinguer deux niveaux : des effets globaux décrivant la progression de la maladie pour l'ensemble de la population, et une description individualisée de la maladie permettant une modélisation spécifique à chaque sujet. Les objectifs des modèles de progression des maladies sont nombreux.

**Diagnostic et détection:** Un espoir est que, en exploitant la quantité considérable de données collectées à partir des cohortes longitudinales, les modèles soient capables de détecter de petits signes de l'apparition de la maladie avant le diagnostic. D'autre part, ils peuvent également être utilisés pour déterminer numériquement les seuils pour le diagnostic de la maladie. De tels modèles pourraient être utilisés pour aider le clinicien dans le diagnostic, et de manière plus réaliste, ils pourraient être mis en place dans les cabinets des médecins généralistes pour suggérer quand orienter le patient vers une consultation chez un neurologue expert.

**Pronostic :** Les modèles de progression des maladies fournissent souvent des prédictions individuelles dans le temps, ce qui permet de prévoir l'évolution des individus dans le futur. Cette tâche est souvent utilisée comme motivation pour la conception d'un modèle de progression des maladies, car la prédiction est plus familière à la communauté de l'apprentissage automatique. Cela a donné

<footnote n="83"><sup>83</sup>C. Murphy, "Olfactory and other sensory impairments in Alzheimer disease," en, *Nature Reviews Neurology*, vol. 15, no. 1, pp. 11–24, Jan. 2019.</footnote>

<footnote n="84"><sup>84</sup>R. B. Postuma, A. E. Lang, J. F. Gagnon, et al., "How does parkinsonism start? Prodromal parkinsonism motor changes in idiopathic REM sleep behaviour disorder," *Brain*, vol. 135, no. 6, pp. 1860–1870, Jun. 2012.</footnote>

<footnote n="85"><sup>85</sup>R. C. Petersen, P. S. Aisen, L. A. Beckett, et al., "Alzheimer's Disease Neuroimaging Initiative (ADNI): Clinical characterization," en, *Neurology*, vol. 74, no. 3, pp. 201–209, Jan. 2010.</footnote>

<footnote n="86"><sup>86</sup>K. Marek, D. Jennings, S. Lasch, et al., "The Parkinson Progression Marker Initiative (PPMI)," en, *Progress in Neurobiology*, Biological Markers for Neurodegenerative Diseases, vol. 95, no. 4, pp. 629–635, Dec. 2011.</footnote>
lieu à la création de défis de prédiction pour la progression des maladies neurodégénératives, tels que le défi TADPOLE<sup>87</sup>. Bien que ce cadre compétitif ait suscité de nombreuses tentatives de modélisation intéressantes, les performances ne devraient pas être le seul objectif d'un modèle.

**Interprétabilité** : Le modèle devrait fournir des informations sur l'évolution typique de la maladie. Être capable d'exposer les paramètres responsables de la spécificité individuelle peut être utilisé pour effectuer une analyse de covariables, par exemple. Une étape importante consiste à concevoir ce que nous appelons une chronologie de la maladie. Il est bien connu que l'âge n'est pas l'axe de référence approprié pour comparer les patients. En effet, certains peuvent commencer la maladie plus tôt, d'autres progressent plus rapidement, ou peut-être les deux en même temps. Deux patients du même âge peuvent ne pas être au même stade de la maladie. D'où la nécessité d'aligner temporellement les patients sur une chronologie commune de la maladie. Une fois cela fait, la variabilité restante dans la présentation de la maladie peut être qualifiée de variabilité spatiale, par opposition à la variabilité temporelle qui tient compte du réalignment des patients dans le temps. Un effort spécifique dans notre travail a été porté sur la distinction entre les effets temporels et spatiaux.

Une fois qu'un modèle de progression des maladies s'est avéré capable d'accomplir de manière fiable l'une des tâches mentionnées précédemment, l'étape suivante consiste à utiliser ce modèle pour atteindre l'objectif ultime de la recherche sur les maladies neurodégénératives : trouver un remède. Un premier espoir qui est rapidement écarté est de trouver une cible potentielle pour un médicament. Cela est très peu probable, car les modèles sont basés sur des données qui ont été collectées pour les besoins spécifiques de la cohorte. Les données sont donc enregistrées en raison de leur lien déjà connu, ou hypothétique, avec le processus neurodégénératif. Une deuxième utilisation potentielle des modèles de progression des maladies est d'identifier de manière plus précise quand, comment et à qui administrer un médicament. Cela est particulièrement utile à l'ère de la médecine personnalisée. Un dernier contexte d'application est l'amélioration des essais cliniques. Récemment, la Food and Drug Administration des États-Unis a approuvé l'utilisation de PROCOVA<sup>88</sup>, une méthode où l'utilisation d'un score de progression de la maladie pronostique dans une analyse de covariance permet un essai plus puissant, ou de manière équivalente, permet une taille d'échantillon plus réduite.

<footnote n="87"><sup>87</sup>R. V. Marinescu, N. P. Oxtoby, A. L. Young, et al., "TADPOLE Challenge: Accurate Alzheimer's disease prediction through crowdsourced forecasting of future data," eng, *PRedictive Intelligence in MEdicine. PRIME (Workshop)*, vol. 11843, pp. 1-10, Oct. 2019.</footnote>

<footnote n="88"><sup>88</sup>D. Bertolini, K. Arnemann, D. Hall, et al., "Machine Learning Enables Smaller ALS Clinical Trials (P1-13.003)," en, *Neurology*, vol. 98, no. 18 Supplement, May 2022.</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 2: Représentation schématique des modèles de progression. Des données de diverses natures répétées au cours de la vie du patient sont collectées dans des bases de données longitudinales qui servent à entraîner des modèles de progression.</figcaption></figure>
Modèles de progression de la maladie

Découlant du modèle hypothétique de<sup>89</sup>, des modèles de progression des maladies ont émergé pour valider l'hypothèse et fournir un modèle basé sur les données en tant que preuve. Une classe spécifique de modèles se concentre sur l'ordre de la progression des anomalies des biomarqueurs. Ces modèles sont appelés modèles basés sur les événements (EBM), et remontent à<sup>90</sup>. Leur application aux maladies neurodégénératives utilise des données transversales pour déterminer l'ordre des événements<sup>91</sup>. Lorsque les marqueurs sont continus, les événements sont définis en fixant des seuils d'anomalie, ce qui permet d'inclure tout type de données dans les modèles. Les EBM ont été étendus pour tenir compte de l'hétérogénéité dans l'algorithme d'inférence des sous-types et des stades (SUSTAIN)<sup>92</sup>. Cependant, la séquence des événements ne fournit que des stades de la maladie, et la chronologie de la maladie n'est pas explicitement modélisée.

De nombreux modèles se sont concentrés sur l'alignement temporel en tant que pierre angulaire pour estimer la progression de la maladie. Un travail précoce a proposé une correspondance linéaire pour extraire la chronologie de la maladie à partir des observations individuelles<sup>93</sup>. Cela a ensuite été amélioré par l'introduction de déformations temporelles<sup>94</sup>. La construction de la chronologie commune de la maladie en faisant correspondre l'âge des patients sur un même axe, par le biais de ce qu'on appelle une reparamétrisation temporelle, considère que les observations des individus sont des aperçus d'un schéma de progression commun global, comme dans<sup>95</sup>. Dans cette étude, la trajectoire moyenne a été estimée sans aucune hypothèse de forme préalable, en utilisant uniquement les aperçus individuels réalignés. Cette chronologie commune de la maladie peut être considérée comme une version continue de la mise en scène de la maladie, un score de progression de la maladie<sup>96</sup>.

La principale classe de modèles longitudinaux utilisés dans la modélisation de la progression des maladies a traditionnellement été la classe des modèles à effets mixtes<sup>97,98</sup>. Ce type de modèles a une

<footnote n="89"><sup>89</sup>C. R. Jack, D. S. Knopman, W. J. Jagust, et al., "Hypothetical model of dynamic biomarkers of the Alzheimer's pathological cascade," English, *The Lancet Neurology*, vol. 9, no. 1, pp. 119-128, Jan. 2010.</footnote>

<footnote n="90"><sup>90</sup>L. A. Beckett, "Maximum Likelihood Estimation in Mallows's Model Using Partially Ranked Data," en, in *Probability Models and Statistical Analyses for Ranking Data*, M. A. Fligner and J. S. Verducci, Eds., ser. Lecture Notes in Statistics, New York, NY: Springer, 1993, pp. 92-107.</footnote>

<footnote n="91"><sup>91</sup>H. M. Fonteijn, M. Modat, M. J. Clarkson, et al., "An event-based model for disease progression and its application in familial Alzheimer's disease and Huntington's disease," eng, *NeuroImage*, vol. 60, no. 3, pp. 1880-1889, Apr. 2012.</footnote>

<footnote n="92"><sup>92</sup>A. L. Young, R. V. Marinescu, N. P. Oxtoby, et al., "Uncovering the heterogeneity and temporal complexity of neurodegenerative diseases with Subtype and Stage Inference," en, *Nature Communications*, vol. 9, no. 1, pp. 1-16, Oct. 2018.</footnote>

<footnote n="93"><sup>93</sup>E. Yang, M. Farnum, V. Lobanov, et al., "Quantifying the pathophysiological timeline of Alzheimer's disease," eng, *Journal of Alzheimer's disease: JAD*, vol. 26, no. 4, pp. 745-753, 2011.</footnote>

<footnote n="94"><sup>94</sup>S. Durrleman, X. Pennec, A. Trouvé, et al., "Toward a Comprehensive Framework for the Spatiotemporal Statistical Analysis of Longitudinal Shape Data," en, *International Journal of Computer Vision*, vol. 103, no. 1, pp. 22-59, May 2013.</footnote>

<footnote n="95"><sup>95</sup>M. C. Donohue, H. Jacqmin-Gadda, M. Le Goff, et al., "Estimating long-term multivariate progression from short-term data," *Alzheimer's & dementia: the journal of the Alzheimer's Association*, vol. 10, no. 0, S400-S410, Oct. 2014.</footnote>

<footnote n="96"><sup>96</sup>V. Kmetzsch, E. Becker, D. Saracino, et al., "Disease Progression Score Estimation From Multimodal Imaging and MicroRNA Data Using Supervised Variational Autoencoders," *IEEE Journal of Biomedical and Health Informatics*, vol. 26, no. 12, pp. 6024-6035, Dec. 2022.</footnote>

<footnote n="97"><sup>97</sup>N. M. Laird and J. H. Ware, "Random-effects models for longitudinal data," eng, *Biometrics*, vol. 38, no. 4, pp. 963-974, Dec. 1982.</footnote>

<footnote n="98"><sup>98</sup>M. J. Lindstrom and D. M. Bates, "Newton-Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data," *Journal of the American Statistical Association*, vol. 83, no. 404, pp. 1014-1022, Dec.</footnote>
structure hiérarchique permettant de différencier les caractéristiques partagées par la population,
ou les effets fixes, de la variabilité individuelle, ou les effets aléatoires. En raison de la nature
non linéaire de la progression de la plupart des biomarqueurs, les modèles sont à effets mixtes
non linéaires<fnref n="99" />. Une famille de courbes de référence est choisie pour la progression, par exemple
des courbes logistiques, et avec une reparamétrisation temporelle affine, les individus sont mis
en correspondance avec la trajectoire moyenne, comme dans<fnref n="100" />. D'autres choix de courbes de
référence sont possibles, comme celles proposées par<fnref n="101" />. Dans le package R pour les modèles à classes
latentes à effets mixtes (LCMM)<fnref n="102" />, le modèle non linéaire peut être combiné avec l'estimation
conjointe des sous-types. Les modèles à effets mixtes non linéaires décrivant la progression des
marqueurs continus peuvent être couplés à la modélisation des événements, dans ce qu'on appelle
des modèles conjoints<fnref n="103" />. Une approche intéressante utilisant des processus gaussiens (GP) pour
tenir compte des effets aléatoires du modèle à effets mixtes<fnref n="104" />,<fnref n="105" />. L'approche Disease Course
Mapping (DCM) utilisée dans<fnref n="106" /> repose sur un modèle géométrique supposant que les observations
se trouvent sur une variété riemannienne. Elle s'inspire d'un modèle géométrique antérieur<fnref n="107" />, qui
était computationnellement inabordable car le bruit était défini dans la variété.

Ces modèles non linéaires ont également été étendus à des modalités de haute dimension telles
que l'imagerie. En prolongement de<fnref n="108" />, les auteurs de<fnref n="109" /> ont appliqué le modèle aux images dans
ADNI. DIVE<fnref n="110" /> est une approche basée sur les voxels pour modéliser l'évolution du cerveau dans
les maladies neurodégénératives. Le modèle GP a été étendu à la modélisation de formes dans<fnref n="111" />.

<footnote n="1988">1988.</footnote>

<footnote n="99">⁹⁹J. Serroyen, G. Molenberghs, G. Verbeke, et al., "Nonlinear Models for Longitudinal Data," *The American Statistician*, vol. 63, no. 4, pp. 378-388, Nov. 2009.</footnote>

<footnote n="100">¹⁰⁰B. M. Jedynak, A. Lang, B. Liu, et al., "A Computational Neurodegenerative Disease Progression Score: Method and Results with the Alzheimer's Disease Neuroimaging Initiative Cohort," *NeuroImage*, vol. 63, no. 3, pp. 1478-1486, Nov. 2012.</footnote>

<footnote n="101">¹⁰¹L. L. Raket, "Statistical Disease Progression Modeling in Alzheimer Disease," *Frontiers in Big Data*, vol. 3, 2020.</footnote>

<footnote n="102">¹⁰²C. Proust-Lima, V. Philipps, and B. Liquet, "Estimation of Extended Mixed Models Using Latent Classes and Latent Processes: The R Package lccmm," en, *Journal of Statistical Software*, vol. 78, pp. 1–56, Jun. 2017.</footnote>

<footnote n="103">¹⁰³B. He and S. Luo, "Joint modeling of multivariate longitudinal measurements and survival data with applications to Parkinson's disease," *eng, Statistical Methods in Medical Research*, vol. 25, no. 4, pp. 1346–1358, Aug. 2016.</footnote>

<footnote n="104">¹⁰⁴M. Lorenzi, X. Pennec, G. B. Frisoni, et al., "Disentangling normal aging from Alzheimer's disease in structural magnetic resonance images," en, *Neurobiology of Aging*, Novel Imaging Biomarkers for Alzheimer's Disease and Related Disorders (NIBAD), vol. 36, S42–S52, Jan. 2015.</footnote>

<footnote n="105">¹⁰⁵M. Lorenzi, M. Filippone, G. B. Frisoni, et al., "Probabilistic disease progression modeling to characterize diagnostic uncertainty: Application to staging and prediction in Alzheimer's disease," en, *NeuroImage*, Mapping diseased brains, vol. 190, pp. 56–68, Apr. 2019.</footnote>

<footnote n="106">¹⁰⁶J.-B. Schiratti, S. Allassonnière, A. Routier, et al., "A Mixed-Effects Model with Time Reparametrization for Longitudinal Univariate Manifold-Valued Data," en, in *Information Processing in Medical Imaging*, S. Ourselin, D. C. Alexander, C.-F. Westin, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2015, pp. 564–575.</footnote>

<footnote n="107">¹⁰⁷T. Fletcher, "Geodesic Regression on Riemannian Manifolds," en, *Sep. 2011*, p. 75.</footnote>

<footnote n="108">¹⁰⁸B. M. Jedynak, A. Lang, B. Liu, et al., "A Computational Neurodegenerative Disease Progression Score: Method and Results with the Alzheimer's Disease Neuroimaging Initiative Cohort," *NeuroImage*, vol. 63, no. 3, pp. 1478-1486, Nov. 2012.</footnote>

<footnote n="109">¹⁰⁹M. Bilgel, J. L. Prince, D. F. Wong, et al., "A multivariate nonlinear mixed effects model for longitudinal image analysis: Application to amyloid imaging," en, *NeuroImage*, vol. 134, pp. 658–670, Jul. 2016.</footnote>

<footnote n="110">¹¹⁰R. V. Marinescu, A. Eshaghi, M. Lorenzi, et al., "DIVE: A spatiotemporal progression model of brain pathology in neurodegenerative disorders," en, *NeuroImage*, vol. 192, pp. 166–177, May 2019.</footnote>

<footnote n="111">¹¹¹C. Abi Nader, N. Ayache, P. Robert, et al., "Monotonic Gaussian Process for spatio-temporal disease progression modeling in brain imaging data," eng, *NeuroImage*, vol. 205, p. 116 266, Jan. 2020.</footnote>
Dans<fnref n="112" />, le modèle DCM a été appliqué à l'imagerie cérébrale, où le cerveau est modélisé comme un maillage intégré dans une variété de formes.

La modélisation de la progression dynamique des biomarqueurs a également été réalisée à l'aide de modèles dynamiques. Ils sont basés sur les équations différentielles ordinaires<fnref n="113" />,<fnref n="114" />,<fnref n="115" />, où le champ des EDO est estimé à partir de l'ensemble de la cohorte et est semblable aux effets fixes des modèles mixtes, tandis que les conditions initiales de la trajectoire d'un patient individuel tiennent compte des effets aléatoires. Les EDO neurales<fnref n="116" /> utilise un réseau d'apprentissage profond pour apprendre le champ. D'autres approches d'apprentissage profond utilisant des réseaux neuronaux récurrents ont également été appliquées à la progression des maladies<fnref n="117" />. Ces méthodes d'apprentissage profond sont souvent moins interprétables. C'est tout le contraire du modèle causal multifactoriel proposé dans<fnref n="118" />. Ce modèle causal permet de réaliser des interventions, ce qui peut être exploité pour émuler une intervention thérapeutique.

Toutes ces méthodes statistiques ont en commun la capacité de décrire la progression de la maladie au niveau de la population. La plupart d'entre elles peuvent également calculer une trajectoire individualisée pour chaque patient. Les principales différences résident dans leur interprétabilité. Les modèles avec une reparamétrisation temporelle explicite fournissent une séparation de la variabilité temporelle et spatiale, ce qui est un atout pour nous. Les modèles semi-paramétriques ou non paramétriques, tels que les modèles dynamiques basés sur les EDO ou les méthodes basées sur l'apprentissage profond, manquent souvent de côté interprétable. Notre travail se concentre sur l'approche géométrique du DCM, qui offre à la fois le désentrelacement spatio-temporel et des paramètres significatifs.

## Maladie de Parkinson

La majeure partie du travail de cette thèse a été développée pour des applications aux données sur la maladie de Parkinson, nous fournirons donc un bref aperçu de cette maladie, de ses causes hypothétiques, de ses manifestations et de ses traitements.

La maladie de Parkinson a été décrite pour la première fois par James Parkinson en 1817. Cette maladie neurodégénérative est la deuxième plus répandue après la maladie d'Alzheimer. Sa prévalence augmente considérablement avec l'âge, passant d'une quasi-nullité avant l'âge de 50 ans à près de 2% après 80 ans<fnref n="119" />. Les hommes ont deux fois plus de chances de développer la maladie

<footnote n="112">¹¹²I. Koval, "Learning Multimodal Digital Models of Disease Progression from Longitudinal Data: Methods & Algorithms for the Description, Prediction and Simulation of Alzheimer's Disease Progression," en, Ph.D. dissertation, Institut Polytechnique de Paris, Jan. 2020.</footnote>

<footnote n="113">¹¹³B. O. Taddé, H. Jacqmin-Gadda, J.-F. Dartigues, et al., "Dynamic modeling of multivariate dimensions and their temporal relationships using latent processes: Application to Alzheimer's disease," en, *Biometrics*, vol. 76, no. 3, pp. 886-899, 2020.</footnote>

<footnote n="114">¹¹⁴N. P. Oxtoby, A. L. Young, D. M. Cash, et al., "Data-driven models of dominantly-inherited Alzheimer's disease progression," *Brain*, vol. 141, no. 5, pp. 1529-1544. May 2018.</footnote>

<footnote n="115">¹¹⁵K. Lahouel, M. Wells, V. Rielly, et al., *Learning nonparametric ordinary differential equations from noisy data*, Feb. 2023.</footnote>

<footnote n="116">¹¹⁶R. T. Q. Chen, Y. Rubanova, J. Bettencourt, et al., *Neural Ordinary Differential Equations*, Dec. 2019.</footnote>

<footnote n="117">¹¹⁷M. Nguyen, T. He, L. An, et al., "Predicting Alzheimer's disease progression using deep recurrent neural networks," en, *NeuroImage*, vol. 222, p. 117-203, Nov. 2020.</footnote>

<footnote n="118">¹¹⁸Y. Iturria-Medina, F. M. Carbonell, R. C. Sotero, et al., "Multifactorial causal model of brain (dis)organization and therapeutic intervention: Application to Alzheimer's disease," en, *NeuroImage*, vol. 152, pp. 60-77, May 2017.</footnote>

<footnote n="119">¹¹⁹A. Elbaz, L. Carcaillon, S. Kab, et al., "Epidemiology of Parkinson's disease," en, *Revue Neurologique*, Neuroepidemiology, vol. 172, no. 1, pp. 14-26, Jan. 2016.</footnote>
que les femmes.

La maladie de Parkinson est caractérisée par la perte de neurones dopaminergiques, qui com-
mence dans la substance noire, une région du tronc cérébral contenant presque exclusivement des
neurones dopaminergiques. La dopamine est un neurotransmetteur impliqué dans le système de
récompense, mais elle est également responsable de l'activation du système moteur. La diminu-
tion des niveaux de dopamine due à la neurodégénérescence entraîne tous les symptômes moteurs
de la maladie de Parkinson que nous connaissons : la bradykinésie (lenteur des mouvements),
la rigidité musculaire, les tremblements, qui sont les principaux symptômes<fnref n="120" />, mais de nombreux
autres symptômes sont également observés, notamment les chutes, la lenteur de la parole... La perte
de neurones dopaminergiques se propage dans le temps à d'autres régions du cerveau, entraînant
des symptômes non moteurs, tels que l'apathie, les hallucinations ou la dépression. Cependant,
l'apparition des symptômes est un indicateur très tardif de la mort neuronale. Lorsque la maladie
de Parkinson est diagnostiquée, déjà 80% des neurones dopaminergiques de la substance noire ont
déjà été perdus. On pense que la neurodégénérescence est due à l'accumulation de la protéine α-
synucléine, qui est actuellement la principale cible du développement de médicaments<fnref n="121" />. La perte
de neurones dopaminergiques peut être observée grâce à la scintigraphie SPECT (tomographie par
émission de photons uniques) du transporteur de dopamine.

Comme mentionné dans les défis liés aux maladies neurodégénératives, la définition de la maladie
de Parkinson est complexe. Elle est actuellement clinique, à la discrétion du neurologue, mais il
n'y a pas de ligne directrice claire pour le diagnostic. Les symptômes sont similaires à ceux de la
maladie à corps de Lewy, ce qui conduit parfois à une révision du diagnostic précoce. Cependant,
nous sommes encore loin d'une définition biologique de la maladie. L'hétérogénéité de la maladie
joue un rôle important, à partir du début des symptômes<fnref n="122" />. Une première distinction peut être
faite en fonction de l'emplacement des premiers signes moteurs. 60% des patients présentent d'abord
une atteinte d'un côté du corps, les symptômes se propageant lentement vers le reste du corps au fil
du temps. Des tentatives de sous-typage ont été faites pour distinguer les sous-types prodromiques
précoces, basés sur les comportements liés au sommeil<fnref n="123" />. L'une des classifications les plus courantes
basées sur les symptômes moteurs est le système dominant le tremblement (TD) par rapport au
trouble de l'équilibre postural et de la marche (PIGD) proposé dans<fnref n="124" />. Cependant, ces sous-types
ne regroupent pas tous les patients, et ils ne sont pas stables dans le temps<fnref n="125" />,<fnref n="126" />. Certains sous-
types particuliers peuvent être dérivés de variants génétiques GBA et LRKK2, qui représentent
environ 5% de l'ensemble des cas de la maladie de Parkinson, les autres étant idiopathiques.

La maladie de Parkinson a bénéficié de la découverte d'un traitement au début des années 1950,
appelé lévodopa. Il s'agit d'un précurseur de la dopamine, et lorsqu'il est pris régulièrement, il com-

<footnote n="120">¹²⁰J. Jankovic, M. McDermott, J. Carter, et al., "Variable expression of Parkinson's disease: A base-line analysis of the DATATOP cohort. The Parkinson Study Group," eng, *Neurology*, vol. 40, no. 10, pp. 1529-1534, Oct. 1990.</footnote>

<footnote n="121">¹²¹C. R. Fields, N. Bengoa-Vergniory, and R. Wade-Martins, "Targeting Alpha-Synuclein as a Therapy for Parkinson's Disease," *Frontiers in Molecular Neuroscience*, vol. 12, 2019.</footnote>

<footnote n="122">¹²²W. J. Zetusky, J. Jankovic, and F. J. Pirozzolo, "The heterogeneity of Parkinson's disease: Clinical and prognostic implications," eng, *Neurology*, vol. 35, no. 4, pp. 522-526, Apr. 1985.</footnote>

<footnote n="123">¹²³R. B. Postuma, A. E. Lang, J. F. Gagnon, et al., "How does parkinsonism start? Prodromal parkinsonism motor changes in idiopathic REM sleep behaviour disorder," *Brain*, vol. 135, no. 6, pp. 1860-1870, Jun. 2012.</footnote>

<footnote n="124">¹²⁴J. Jankovic, M. McDermott, J. Carter, et al., "Variable expression of Parkinson's disease: A base-line analysis of the DATATOP cohort. The Parkinson Study Group," eng, *Neurology*, vol. 40, no. 10, pp. 1529-1534, Oct. 1990.</footnote>

<footnote n="125">¹²⁵T. Simuni, C. Caspell-Garcia, C. Coffey, et al., "How stable are Parkinson's disease subtypes in de novo patients: Analysis of the PPMI cohort?" en, *Parkinsonism & Related Disorders*, vol. 28, pp. 62-67, Jul. 2016.</footnote>

<footnote n="126">¹²⁶R. S. Eisinger, C. W. Hess, D. Martinez-Ramirez, et al., "Motor subtype changes in early Parkinson's disease," en, *Parkinsonism & Related Disorders*, vol. 43, pp. 67-72, Oct. 2017.</footnote>
pense le manque de production de dopamine dans le cerveau. D'autres médicaments à base d'apport dopaminergique ont depuis été développés, tels que les agonistes dopaminergiques ou les inhibiteurs de recapture, notamment IMAOB et ICOMT. Ils reposent tous sur le même mécanisme, avec des effets secondaires différents. Par exemple, on sait que la lévodopa provoque des dyskinésies, c'est-à-dire des mouvements spontanés, et a également été liée à des blocages dans certaines études<fnref n="127" />. Les agonistes dopaminergiques ont plus de chances de provoquer des troubles du contrôle des impulsions (TCI), qui entraînent des comportements addictifs tels que le jeu. Le traitement des patients atteints de la maladie de Parkinson consiste aujourd'hui en un cocktail de plusieurs médicaments dopaminergiques, et les experts débattent toujours de la meilleure politique de traitement<fnref n="128" />. De nombreux facteurs peuvent influencer la réponse d'un patient au traitement et sa probabilité de développer des effets secondaires négatifs.

Le traitement dopaminergique est purement symptomatique, il n'a montré aucune amélioration à long terme<fnref n="129" />. Son efficacité diminue à mesure que la maladie progresse vers des stades plus avancés, entraînant des fluctuations entre les états "ON" et "OFF". Les états "ON" correspondent à l'état de traitement fonctionnant normalement, où les symptômes moteurs sont atténués. Les états "OFF", en revanche, surviennent lorsque l'effet du traitement s'estompe. Aux stades avancés de la maladie de Parkinson, la demi-vie du médicament dopaminergique devient très courte en raison de la quasi-absence de production par le cerveau.

Un autre traitement appelé stimulation cérébrale profonde consiste à implanter une électrode dans l'amygdale pour stimuler directement le cerveau. Ce traitement s'est révélé être une alternative puissante à la thérapie dopaminergique<fnref n="130" />, mais il reste assez rare en raison de son caractère invasif. Le traitement dopaminergique a longtemps éclipsé la recherche d'un médicament modifiant la maladie, ce qui fait que la recherche dans ce domaine est bien inférieure à celle de la maladie d'Alzheimer.

Les cohortes longitudinales utilisées dans ce travail de thèse comprennent PPMI<fnref n="131" /> et une co-horte française du consortium NS-PARK. Les données des cohortes sur la maladie de Parkinson sont principalement composées d'évaluations cognitives et motrices. Historiquement, l'échelle de Hoehn et Yahr<fnref n="132" /> (de 0 à 5) était utilisée pour décrire la maladie selon des stades successifs. Aujourd'hui, l'échelle la plus fréquemment utilisée pour évaluer la progression de la maladie est l'échelle unifiée de la maladie de Parkinson du Mouvement Disorder Society (MDS-UPDRS)<fnref n="133" />, qui regroupe de nombreux éléments évaluant différents aspects de la maladie, notamment un examen moteur et une évaluation de l'impact de la maladie sur la vie quotidienne du patient.

<footnote n="127">¹²⁷M. Gilat, N. D'Cruz, P. Ginis, et al., "Freezing of gait and levodopa," English, *The Lancet Neurology*, vol. 20, no. 7, pp. 505–506, Jul. 2021.</footnote>

<footnote n="128">¹²⁸R. M. A. de Bie, C. E. Clarke, A. J. Espay, et al., "Initiation of pharmacological therapy in Parkinson's disease: When, why, and how," eng, *The Lancet. Neurology*, vol. 19, no. 5, pp. 452–461, May 2020.</footnote>

<footnote n="129">¹²⁹R. Cilia, E. Cereda, A. Akpalu, et al., "Natural history of motor symptoms in Parkinson's disease and the long-duration response to levodopa," *Brain*, vol. 143, no. 8, pp. 2490–2501, Aug. 2020.</footnote>

<footnote n="130">¹³⁰A. L. Benabid, "Deep brain stimulation for Parkinson's disease," en, *Current Opinion in Neurobiology*, vol. 13, no. 6, pp. 696–706, Dec. 2003.</footnote>

<footnote n="131">¹³¹K. Marek, D. Jennings, S. Lasch, et al., "The Parkinson Progression Marker Initiative (PPMI)," en, *Progress in Neurobiology*, Biological Markers for Neurodegenerative Diseases, vol. 95, no. 4, pp. 629–635, Dec. 2011.</footnote>

<footnote n="132">¹³²M. M. Hoehn and M. D. Yahr, "Parkinsonism: Onset, progression, and mortality," en, *Neurology*, vol. 17, no. 5, pp. 427–427, May 1967.</footnote>

<footnote n="133">¹³³C. G. Goetz, B. C. Tilley, S. R. Shaftman, et al., "Movement Disorder Society-sponsored revision of the Unified Parkinson's Disease Rating Scale (MDS-UPDRS): Scale presentation and clinimetric testing results: MDS-UPDRS: Clinimetric Assessment," en, *Movement Disorders*, vol. 23, no. 15, pp. 2129–2170, Nov. 2008.</footnote>
## Buts de la thèse

Le travail de cette thèse s'est appuyé sur le modèle de Disease Course Mapping (DCM) introduit dans<fnref n="134" />. Nous avons décrit certains des principaux défis posés par les maladies neurodégénératives et nous allons étendre le modèle DCM pour relever ces défis. Le modèle DCM a déjà été conçu pour tenir compte de l'alignement temporel et de la séparation spatio-temporelle. Cependant, sa conception est loin d'être parfaite. Tout d'abord, il est limité aux données continues. De plus, le modèle DCM ne tient pas compte de l'hétérogénéité des sous-types de maladies. Enfin, il n'est pas adapté aux applications dans le développement de médicaments. Notre bref aperçu de la maladie de Parkinson met en évidence ces limitations. Les symptômes, la définition de stades de la maladie et les échelles de notation sont omniprésents, et toutes les variables sont non-continues. Les sous-types sont essentiels pour modéliser la variabilité de la maladie de Parkinson. Enfin, toutes les données des cohortes sur la maladie de Parkinson proviennent de patients sous traitement dopaminergique, ce qui nécessite un modèle spécifique pour cela. Nous fournirons donc les extensions suivantes du modèle DCM :

* L'applicabilité aux données discrètes, y compris les observations binaires et ordinales.

* La structure de mélange pour découvrir les sous-types et mieux expliquer l'hétérogénéité de la maladie.

* Des modèles de traitement liés au modèle DCM pour prendre en compte l'effet du traitement et ouvrir la voie au développement de médicaments modificateurs de la maladie.

Toutes nos contributions méthodologiques sont motivées par des besoins pratiques rencontrés dans des données réelles. Elles sont donc accompagnées d'applications à des cohortes neurodégénératives. Nous avons utilisé l'interprétabilité du modèle pour fournir des informations cliniques dans l'analyse de nos résultats expérimentaux.

## Plan du manuscrit

Le premier chapitre introduit le modèle Disease Course Mapping (DCM) et discute de ses limites. La deuxième partie du chapitre présente Geodesic Bending, une extension du modèle DCM qui introduit plus de flexibilité dans le modèle en apprenant la forme de la trajectoire typique de la maladie.

Le deuxième chapitre présente notre extension du modèle DCM aux données binaires et ordinaires. Nous fournissons des applications approfondies aux cohortes de la maladie de Parkinson, ainsi qu'une mise en parallèle avec la théorie de la réponse à l'item.

Le troisième chapitre aborde le défi de l'hétérogénéité des maladies neurodégénératives. Notre contribution prend la forme de deux variantes de mélange pour le modèle DCM. Nous mettons en évidence les différences entre les deux avec une application à la maladie d'Alzheimer.

Le quatrième chapitre se concentre sur la modélisation du traitement. Nous présentons deux approches de modélisation. La première consiste en un cadre générique utilisant n'importe quel

<footnote n="134">¹³⁴J.-B. Schiratti, S. Allassonnière, A. Routier, et al., "A Mixed-Effects Model with Time Reparametrization for Longitudinal Univariate Manifold-Valued Data," en, in *Information Processing in Medical Imaging*, S. Ourselin, D. C. Alexander, C.-F. Westin, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2015, pp. 564-575.</footnote>
modèle de progression de la maladie pour estimer l'effet du traitement à partir des seules données des patients, sans nécessiter de données de témoins ni de placebo. La deuxième approche s'appelle le modèle géodésique par morceaux et est destinée aux médicaments modificateurs de la maladie, fournissant un cadre pour apprendre les changements dans la trajectoire de la maladie dus au traitement.

Le cinquième chapitre résume toutes nos contributions à l’algorithme d’estimation dans toutes les extensions précédentes. Il présente l’algorithme MCMC-SAEM utilisé pour l’estimation du DCM. Nous introduisons ensuite nos contributions sous la forme de variantes spécifiques de cet algorithme.
<figure><img src="image_9.png" /></figure>
# Chapter 1

## The Disease Course Mapping model: limits and first extension

This chapter will introduce disease course mapping (DCM), a particular disease modeling technique. This model was initially introduced in the original paper<fnref n="1" />. We refer the reader to the related thesis work that provides an in-depth presentation of the model for further details<fnref n="2" />. The model was later extended and successfully applied to neurodegenerative diseases by subsequent PhD students in their thesis work<fnref n="3" />,<fnref n="4" />,<fnref n="5" />. We will describe the structure of the model and provide an intuition of its parameters and geometric interpretation. Then we will focus on a specific variant of the model focused on logistic curves, as they will be used throughout the entirety of this thesis. All the materials about model parameters estimation is tackled in a dedicated chapter 5, as we will discuss the algorithm and all the contributions that we brought to it in the same section. This choice was made purposefully so that the focus in these first chapters is on the model itself, the intuition behind its geometry, the interpretation of the parameters and the applications.

The second part of this chapter concentrates on an extension of the DCM model called Geodesic Bending. This extension addresses one of the limitations of the DCM model, namely the requirement to select a family of template curves. We demonstrate how this limitation can be overcome by learning the shape of disease progression trajectories.

Our contributions in this chapter are the development of the Geodesic Bending model, in collaboration with Samuel Gruffaz. This lead to the paper "Learning Riemannian metric for disease progression modeling"<fnref n="6" />.

<footnote n="1">¹J.-B. Schiratti, S. Allassonnière, O. Colliot, et al., "A Bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations," *The Journal of Machine Learning Research*, vol. 18, no. 1, pp. 4840-4872, Jan. 2017.</footnote>

<footnote n="2">²J.-B. Schiratti, "Methods and algorithms to learn spatio-temporal changes from longitudinal manifold-valued observations," en, Ph.D. dissertation, Université Paris Saclay (COmUE), Jan. 2017.</footnote>

<footnote n="3">³I. Koval, "Learning Multimodal Digital Models of Disease Progression from Longitudinal Data: Methods & Algorithms for the Description, Prediction and Simulation of Alzheimer's Disease Progression," en, Ph.D. dissertation, Institut Polytechnique de Paris, Jan. 2020.</footnote>

<footnote n="4">⁴A. Bône, "Learning adapted coordinate systems for the statistical analysis of anatomical shapes. Applications to Alzheimer's disease progression modeling," en, Ph.D. dissertation, Sorbonne Université, Jan. 2020.</footnote>

<footnote n="5">⁵R. Couronné, "Progression models for Parkinson's Disease," en, Ph.D. dissertation, Sorbonne Université, Sep. 2021.</footnote>

<footnote n="6">⁶S. Gruffaz, P.-E. Poulet, E. Maheux, et al., "Learning Riemannian metric for disease progression modeling," in</footnote>
Contents of the chapter

<table><tr><td>1.1</td><td>The Disease Course Mapping model</td><td>29</td></tr><tr><td>1.1.1</td><td>Notations and mixed-effect modelling</td><td>29</td></tr><tr><td>1.1.2</td><td>Generic disease course mapping framework</td><td>32</td></tr><tr><td>1.1.3</td><td>Bayesian model</td><td>36</td></tr><tr><td>1.1.4</td><td>Limitations</td><td>38</td></tr><tr><td>1.2</td><td>Geodesic Bending</td><td>38</td></tr><tr><td>1.2.1</td><td>Introduction and related work</td><td>39</td></tr><tr><td>1.2.2</td><td>Method</td><td>40</td></tr><tr><td>1.2.3</td><td>Experiments</td><td>43</td></tr><tr><td>1.2.4</td><td>Discussion on Geodesic Bending</td><td>50</td></tr><tr><td>1.3</td><td>Conclusion</td><td>50</td></tr></table>

<footnote>Advances in Neural Information Processing Systems, M. Ranzato, A. Beygelzimer, Y. Dauphin, et al., Eds., vol. 34, Curran Associates, Inc., 2021, pp. 23780–23792.</footnote>
# 1.1 The Disease Course Mapping model

## 1.1.1 Notations and mixed-effect modelling

### Notations

We provide here the introduction of our standard notations for a longitudinal dataset $\mathcal{D}$ and the model trajectories that will be used throughout this work. Our longitudinal dataset is composed of several measurements ($y_{ijk}$)$_{1\le i\le N, 1\le j\le N, 1\le k\le d}$ where $d$ features were observed on the patient $i$ at $N_i$ visits. We often aggregate these features as the coordinates of a single vector in a multidimensional space $\mathbf{y}_{ij}$. $t_{ij}$ is the age of the patient $i$ at the $j$th visit. Some observations can be missing, as the probabilistic framework of the model will allows us to ignore those. The number of visits $N_i$ per patient is variable. Usually observations $y_{ijk}$ are continuous values. We will use the notation $\eta$ specifically for the model's predicted trajectory in a continuous space. With a longitudinal dataset, $\eta$ will always be a function of time $t$ such that the observations follow a distribution $p$ depending on this model $\mathbf{y}_{ij} \sim p(\eta(t_{ij}; \theta))$, where $\theta$ generally refer to model parameters. $p$ will be called noise model. We call trajectory the curve : $t \mapsto \eta(t; \theta)$. The model can use additional covariates $\mathbf{X}_{ij}$ for individual $i$, which may depend on the time of the visit $t_{ij}$.

Let us introduce a simple model to use these notations.

### Linear mixed effect model

Trying to model the observations as a random gaussian noise on top of a simple function of time $y_{ijk} = \eta(t_{ij}) + \epsilon_{ij}$ with $\epsilon_{ij} \sim \mathcal{N}(0, \sigma^2)$ is too simple since the independence assumption of the noise model between all observations does not hold. Indeed the successive observations of a same individual $i$ are correlated, and this structure in the observations should be taken into account. This is where mixed-effect models come into play.

Mixed-effect modelling<sup>7</sup> is the gold standard method to deal with longitudinal data. In this framework, the model $\eta$ has two types of parameters $\theta = (\mathbf{z}_{pop}, (\mathbf{z}_i)_{i\in[1,N]})$. The first set of parameters $\mathbf{z}_{pop}$ represent common effects across the population while $\mathbf{z}_i$ will relate to model parameters that are specific to individual $i$. $\mathbf{z}_{pop}$ are called fixed effects, and we will equally refer to them as *population parameters*. $\mathbf{z}_i$ are called random effects, or alternatively *individual parameters*. Taking the previous example as a starting point, we would have the model $y_{ijk} = \eta(t_{ij}; \mathbf{z}_{pop}, \mathbf{z}_i) + \epsilon_{ij}$ with $\epsilon_{ij} \sim \mathcal{N}(0, \sigma^2)$ and $\eta$ parametrized by both fixed and random effects.

From a conceptual point of view, this type of modelling allows to distinguish between two levels. The population level is defined only by fixed effects, and the longitudinal trajectory where only fixed effects are kept, i.e. $\mathbf{z}_i = 0$, should be representative of the average trend over all the individuals. We often call $\gamma_0 : t \mapsto \eta(t; \mathbf{z}_{pop}, 0)$ the average trajectory. For neurodegenerative diseases, the average trajectory $\gamma_0$ should reflect the knowledge of the typical disease progression. These fixed effects therefore should be more important in terms of total variance explained compared to the random effects. The individual parameters on the other hand are here to adapt the average trajectory to each individual trajectory. These random effects act as slight perturbations around the average trajectory and allow for a better description at an individual level. Again, in the context of neurodegenerative diseases, this is crucial for personalized medicine by adapting the treatment to an individual specific disease progression.

<footnote n="7"><sup>7</sup>R. A. Fisher, "XV.—The Correlation between Relatives on the Supposition of Mendelian Inheritance.," en, *Earth and Environmental Science Transactions of The Royal Society of Edinburgh*, vol. 52, no. 2, pp. 399–433, Jan. 1919.</footnote>
A linear mixed-effect model<fnref n="8" /> is a simple longitudinal mixed-effect model which will be our reference to get the intuition on how the disease course mapping model works. We assume that observations are continuous in the first place, $y_{ij} \in \mathbb{R}^d$, and that the noise is Gaussian. Let us write the model linear trajectory and its noise model as follows:

$$ \eta(t_{ij}; \mathbf{X}_{ij}, \mathbf{z}_{pop}, \mathbf{z}_i) = C \mathbf{X}_{ij} + (\mathbf{v} + \mathbf{v}_i)t_{ij} + (\boldsymbol{\beta} + \boldsymbol{\beta}_i) \quad (1.1.1) $$

$$ \mathbf{y}_{ij} \sim \mathcal{N}(\eta(t_{ij}; \mathbf{z}_{pop}, \mathbf{z}_i), \Sigma) \quad (1.1.2) $$

with the fixed effects (population parameters) $\mathbf{z}_{pop} = (\mathbf{v}, \boldsymbol{\beta}, \mathbf{C})$, the random effects (individual parameters) $\mathbf{z}_i = (\mathbf{v}_i, \boldsymbol{\beta}_i)$ and the noise covariance matrix $\Sigma$ which is often assumed to be diagonal. The random effects are assumed to follow a centered distribution, typically Gaussian $\mathbf{z}_i \sim \mathcal{N}(0, \boldsymbol{\Omega})$ where $\boldsymbol{\Omega}$ is the covariance matrix of the random effects. $\boldsymbol{\Omega}$ can be fixed, or learned with sometimes a specified structure (generally diagonal) in which case $\boldsymbol{\Omega}$ is part of the model parameters $\boldsymbol{\theta}$.

Having in mind the challenges presented in the introduction about spatiotemporal disentangling, we introduce a new formulation of the model above. We keep the same dimension for the parameters, but we separate time-related effects from spatial effects. This means that the random intercept $\boldsymbol{\beta}_i$ is replaced by two new terms $-\tau_i \mathbf{v}$ and $\mathbf{w}_i \in \text{Span}(\mathbf{v})^{\perp}$. We also decompose the random slope in a multiplicative component $\alpha_i$ on $\mathbf{v}$ and a random slope change $\zeta_i$ in $\text{Span}(\mathbf{v})^{\perp}$. The new formula is:

$$ \eta(t_{ij}; \mathbf{X}_{ij}, \mathbf{z}_{pop}, \mathbf{z}_i) = C \mathbf{X}_{ij} + \boldsymbol{\beta} + (\alpha_i t_{ij} - \tau_i)(\mathbf{v} + \boldsymbol{\zeta}_i) + \mathbf{w}_i \quad (1.1.3) $$

We thus have three terms that are related to identifiable effects:

* the covariate effect $\mathbf{X}_{ij}^T \mathbf{C}$ reflects how the trajectory $\eta$ will be influenced by the covariates. Note that the effect for a given covariate $X_{ijk}$ can also be decomposed into a temporal and a spatial component, by decomposing the $k$-th column of the matrix $\mathbf{C}$ on $\text{Span}(\mathbf{v})$ and $\text{Span}(\mathbf{v})^{\perp}$

* the temporal effect $\alpha_i t_{ij} - \tau_i$ such that the random temporal effect is an affine time-warping where we can identify the acceleration factor $\alpha_i$ and the time intercept $\tau_i$. These two parameters do not affect the trajectory in $\mathbb{R}^d$, they only modify the speed of progression and the starting time

* the spatial effects $\zeta_i$ and $\mathbf{w}_i$ which corresponds to changes of the trajectory in $\mathbb{R}^d$, where $\mathbf{w}_i$ is a random shift of the curve parallel to the average trajectory defined by the fixed effects and $\zeta_i$ is a random shift of direction of the trajectory. These correspond to the intercept and slope random effects which are not related to the axis of time

For the rest of this chapter we will focus on the last two terms, assuming that we have no covariates in the model. Modelling with covariates is the purpose of another ongoing work by a colleague<fnref n="9" />. We will therefore simplify the equations hereafter by positing that $\mathbf{X}_{ij} = 0$. We introduce the time-warping function $\psi_i(t) = \alpha_i t - \tau_i$ and the spatial curve function $\gamma_i(t) = \boldsymbol{\beta} + t(\mathbf{v} + \boldsymbol{\zeta}_i) + \mathbf{w}_i$ such that the model becomes $\eta(t; \mathbf{z}_{pop}, \mathbf{z}_i) = \gamma_i(\psi_i(t))$. Notice that $\gamma_0$ corresponds to a spatial curve with $\mathbf{z}_i = 0$.

Figure 1.1 illustrates these two effects.

<footnote n="8">⁸N. M. Laird and J. H. Ware, “Random-effects models for longitudinal data,” eng, Biometrics, vol. 38, no. 4, pp. 963–974, Dec. 1982.</footnote>

<footnote n="9">⁹N. Fournier and S. Durrleman, “Covariate-Aware Longitudinal Modelling for Neurodegenerative Diseases,” en,</footnote>
<figure><img src="image_10.png" /><figcaption>Figure 1.1: Spatial and temporal effects on the trajectory. Left plot: spatial effect, $\gamma_{pop}$ is the mean trajectory defined by population parameters while $\gamma_i$ is the individualized trajectory of individual $i$. Middle and right plot: temporal effect, showing successively the effect of time-shifting with $\tau_i$ and acceleration $\alpha_i$ on $\gamma_i$ to obtain $\eta_i$.</figcaption></figure>

As mentioned previously, trajectories for neurodegenerative markers tend to be non-linear. Therefore we need to introduce a non-linearity in the previous model.

## Non-linear mixed effect model

One option would be to simply apply a non-linear deformation to the linear model above. This is the approach of generalized linear models (GLM)<fnref n="10" /> and one that is often taken in non-linear mixed-effect modelling. For instance the LCMM library uses the linear mixed-effect model as a latent variable which is then passed through a non-linear link function such as a spline or a logit. If we denote $\Lambda$ the previous linear model, the model trajectory would become $\eta(t; \theta_f, z_{pop}, z_i) = f(\Lambda(t; z_{pop}, z_i); \theta_f)$ where the link function $f$ is parametrized by $\theta_f$.

This approach is perfectly sound from a theoretical and a practical point of view. However the use of a non-linear function on top of random effects which are normally distributed in the Euclidean latent space might lead to an awkward distribution in the observation space. Especially when small random effects on the slope in the linear case can lead to large changes in the long run, which can be further amplified by the non-linear deformation, leaving the distribution of random-effects $\zeta_i$ very sensitive to the choice of link function and prone to cause optimization issues. This remark motivated the geometric model developed by<fnref n="11" /> which lead to the disease course mapping model.

<footnote n="10">¹⁰P. McCullagh, *Generalized Linear Models*, en. Routledge, Oct. 2018.</footnote>

<footnote n="11">¹¹J.-B. Schiratti, S. Allassonnière, A. Routier, et al., “A Mixed-Effects Model with Time Reparametrization for Longitudinal Univariate Manifold-Valued Data,” en, in *Information Processing in Medical Imaging*, S. Ourselin, D. C. Alexander, C.-F. Westin, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2015, pp. 564–575.</footnote>
1.1.2 Generic disease course mapping framework

The disease course mapping model (DCM) is a non-linear mixed-effect model based on a geometric model which we describe below. The geometric model is based on Riemannian geometry, and the spatial variability is defined directly in the observation space instead of a latent space. This allows for a better control of the random effects distributions while also being more directly interpretable.

We will not delve into the concrete mathematical notions about the geometric model in this section. This has already been done in previous theses<fnref n="12" />,<fnref n="13" />,<fnref n="14" />,<fnref n="15" />. We prefer to provide an intuition about the model. Some introductory notions of Riemannian geometry are provided in the appendix 15. For more details about Riemannian geometry, we refer the reader to the excellent book<fnref n="16" />. All mathematical formulas for the DCM model used in this thesis are also found in the appendix 15, which provides a concrete summary of the notions presented in this chapter.

Model definition and intuition

The observations $y_{ij}$ are assumed to be points belonging to a Riemannian manifold ($M$, $g$). A Riemannian manifold is a manifold $M$ equipped with a metric $g$. The model considers that repeated observations in time for a subject $i$ follow a trajectory $\gamma_i$ on the manifold. The average trajectory $\gamma_0$ is assumed to be a geodesic in the manifold, which is the equivalent of a straight line in the Euclidean setting. Indeed a geodesic $\gamma$ is the only solution to the no-acceleration differential equation $\nabla_{\gamma}\dot{\gamma} = 0$ with specified initial conditions for $\gamma(t_0)$ and $\dot{\gamma}(t_0)$ (where $\nabla$ is the Levi-Civita connection, which is equivalent to the notion of derivation in a Riemannian manifold).

The DCM model builds individual trajectories $\gamma_i$ around the average trajectory $\gamma_0$. The set of all individual trajectories form a bundle of trajectories in the manifold. The random effects control the divergence of the trajectories in the bundle. In order to avoid trajectories diverging too much from each other in the space, we add a crucial hypothesis: the parallel trajectories. The DCM model stipulates that the random effects on the slopes $\zeta_i$ are negligible and should be removed. Therefore the spatial variability boils down to the shift of the curve $\gamma_0$ parallel to itself on the manifold. The notion of parallels in the DCM has been defined thanks to the Exp-parallelization 10, which extends parallels to Riemannian manifolds. The set of individual trajectories $\gamma_i$ is therefore a bundle of parallel trajectories in a Riemannian sense. However note that the Exp-parallelization does not conserve the dynamic properties of the trajectory $\gamma_0$. Even though $\gamma_0$ is a geodesic, there is no guarantee that parallel trajectories $\gamma_i$ are also geodesics. Figure 1.2 schematizes the process of Exp-parallelization.

This modelling implies that the choice of the Riemannian manifold and metric defines the
possible trajectories. One simple choice is to use $\mathbb{R}^d$ equipped with the euclidean metric, where all
the trajectories will be simple parallel straight lines. In Figure 1.3, one can see another choice of
manifold and metric, leading to trajectories being logistic curves.

<footnote n="12">¹²J.-B. Schiratti, “Methods and algorithms to learn spatio-temporal changes from longitudinal manifold-valued observations,” en, Ph.D. dissertation, Université Paris Saclay (COMUE), Jan. 2017.</footnote>

<footnote n="13">¹³A. Bône, “Learning adapted coordinate systems for the statistical analysis of anatomical shapes. Applications to Alzheimer’s disease progression modeling,” en, Ph.D. dissertation, Sorbonne Université, Jan. 2020.</footnote>

<footnote n="14">¹⁴I. Koval, “Learning Multimodal Digital Models of Disease Progression from Longitudinal Data : Methods & Algorithms for the Description, Prediction and Simulation of Alzheimer’s Disease Progression,” en, Ph.D. dissertation, Institut Polytechnique de Paris, Jan. 2020.</footnote>

<footnote n="15">¹⁵R. Couronné, “Progression models for Parkinson’s Disease,” en, Ph.D. dissertation, Sorbonne Université, Sep. 2021.</footnote>

<footnote n="16">¹⁶M. Do Carmo, *Riemannian Geometry*, en.</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 1.2: Scheme of the Exp-parallelization on a Riemannian manifold. $\gamma_0$ is the reference geodesic and $\gamma_i$ is the trajectory constructed parallely to $\gamma_0$ with the space-shift $\mathbf{w}_i$. $\mathbf{w}_i$ belongs to the tangent plane at $\mathbf{p}$, and we obtain the first point of $\gamma_i$ by following the geodesic starting from $\mathbf{p}$ with derivative $\mathbf{w}_i$ (this operation is called shooting). The other points are obtained by shifting $\mathbf{w}_i$ along the $\gamma_0$ curve using parallel transport. More specifically this is done by shooting from $\gamma_0(t)$ with $P_{(\gamma_0,t_0,t)}(\mathbf{w}_i)$ which is the vector $\mathbf{w}_i$ transported from $t_0$ to $t$ along the $\gamma_0$ curve.</figcaption></figure>

In order to parametrize the distribution of those trajectories together, the model has a mixed-effects hierarchical structure that separates the population fixed parameters and the individual random parameters:

* **Population parameters:** the individual trajectories form a bundle centered around an "average" trajectory. This central trajectory is taken as a geodesic $\gamma_0$ in the manifold $\mathcal{M}$. This geodesic is parametrized by its initial conditions at time $t_0$: the initial point on the manifold $\gamma_0(t_0) = \mathbf{p}$ (thus $\mathbf{p} \in \mathbb{R}^d$) and the initial speed $\dot{\gamma}_0(t_0) = \mathbf{v}$ (of dimension $d$ as well).

* **Individual parameters:** here we separate temporal and spatial parameters

- **Spatial parameters:** each individual trajectory $\gamma_i$ stems from the central geodesic $\gamma_0$ with a small shift on the manifold $\mathbf{w}_i$ called a space-shift, or inter-marker spacing. The individual follows a trajectory parallel to $\gamma_0$ in the sense of the Exp-parallelisation. The space-shift $\mathbf{w}_i$ is a vector in the tangent space at $\mathbf{p}$ (thus of dimension $d$) defined so that shooting with the Riemannian exponential from $\mathbf{p}$ at speed $\mathbf{w}_i$ yields a point of $\gamma_i$, cf Figure 1.2.

- **Temporal parameters:** The individual temporality of the disease can be modelled with a time reparametrization $\psi_i(t)$ so that the observation $\mathbf{y}_{ij}$ is $\gamma_i(\psi_i(t_{ij}))$. The chosen time-warp is $\psi_i(t) = \alpha_i(t - t_0 - \tau_i) + t_0$, where $\tau_i$ is called the time-shift and models an early or late onset, and $\alpha_i = e^{\xi_i}$ is called the acceleration factor.
<figure><img src="image_12.png" /><figcaption>Figure 1.3: **A.** Generic model scheme. The average population trajectory $\gamma_0$ is a geodesic in the manifold. An individual trajectory $\gamma_i$ is defined as an exp-parallellisation of $\gamma_0$ by a space-shift $\mathbf{w}_i$. **p** and **v** are shown on $\gamma_0$, with a space-shift $\mathbf{w}_i$ pointing to the resulting trajectory $\gamma_i$ upon which lie the measurements $(\mathbf{y}_{ij})_j$. **B.** Logistic curves model. A bundle of trajectories corresponding to the logistic curves model in the 2-dimensional manifold is shown. **C.** Personalization process in the logistic curves model. The plots show the successive influence of the individual parameters on the average trajectory (left) until the trajectory is personalized and matches the measurements (right). Note that the x-axis introduces the time and the two dimensions of the measurements are shown as two different colors.</figcaption></figure>

Figure 1.3 provides a visual understanding of the geometric model.

Compared to Figure 1.1, we see in that the spatial variability does not include the random effect on the slope $\zeta_i$ in sub-figure A. Therefore we only have "parallel" trajectories, and the bundle resulting from it is shown in the sub-figure B. Sub-figure C illustrates the temporal effect, which is similar to the one observed in the linear case. The small difference is that the reference of the time axis is centered on $t_0$ rather than the arbitrary 0 in Figure 1.1.

For the spatial parameters, the space-shifts $(\mathbf{w}_i)_i$ have the same dimension as the observations and are a vector in the tangent space at $\mathbf{p}$, as does $\mathbf{v}$. In order to have identifiability with the time delay $\tau_i$ which operates as a shift along the curve, the space-shifts are required to be orthogonal to $\mathbf{v}$, which we enforce by restraining $\mathbf{w}_i$ to $\text{Span}(\mathbf{v})^{\perp}$, where the orthogonality is defined with the
inner product of the Riemannian manifold, which is defined for two vectors **x**, **y** ∈ ℛM_p in the
tangent space of a given point **p** depending on the metric at this point ⟨**x**, **y**⟩_p = g_p(u, v).

For better interpretability the model reduces the dimension of the space-shift parameter. It is
estimated with an independent component analysis (ICA decomposition)<fnref n="17" /> with N_s independent
sources (s_i)₁≤i≤N_s. The formulation unfolds as w_i = As_i where the columns of A are orthogonal
to v. The columns of the matrix A are defined as A_l = Σ_{k=1}^{d-1} β_{lk} B_k a linear combination of an
orthonormal basis (B_k)₁≤k≤d-1 of the orthogonal hyperplane to Span(v). This orthonormal basis is
computed using a Householder method. More details on the formulas are provided in the appendix
.045.

Note that the model is not structurally identifiable: there is indeed an infinity of triplets (t₀, **p**, **v**)
defining the same geodesic γ₀. At t₀ fixed, the parameters (p, v) are unique, and so are the individual
time parameters (τᵢ, ξᵢ). The space related parameters **A** and **s**ᵢ are unique up to a permutation of
the columns of the matrix **A**, parametrized by the (βₗₖ)ₗₖ. Identifiability comes with the statistical
model.

The logistic curves model

The disease course mapping model is a generic framework, and we can apply it to different types
of observations by choosing the adequate manifold. For instance one can model symmetric definite
positive matrices<fnref n="18" /> or meshes<fnref n="19" />. For every choice, one has to compute the metric, the geodesic
formulas and the Exp-parallelization formulas. However in most Riemannian manifolds, the Exp-
parallelization does not yield an explicit formula, thus requiring heavy computations for each indi-
vidual trajectory γᵢ. This limits greatly the possible choices of manifolds. The most useful ones are
the real-valued intervals manifolds with the metric resulting from the push-forward of a function f.
In this case the resulting geodesics are the result of a straight line to which we applied the function
f. For instance we can take the push-forward of the logit, or the exponential function, and the
resulting geodesics would be logistics and exponential curves respectively.

Notice how we end up coming back to the simpler formulation of a non-linear mixed-effect model
initially proposed 1.1.1, which consisted in applying a non-linearity to a linear mixed-effect model.
The result of us using this fancy geometric framework is that the parametrization of the slope **v**
and the space-shifts **w**ᵢ will be adapted to the geometry of the data, with the distribution of the
space-shifts having desirable properties in the manifold.

We will focus on the manifold whose geodesics are logistic curves. The choice of logistic curves also echoes a well-known hypothesis about the progression of biomarkers in Alzheimer's disease<fnref n="20" />, but the use of logistic curves is common in most neurodegenerative diseases. This is very convenient as we will see, since the logistic is also fundamental for discrete data, which we will tackle in the next chapter. All starts with the one-dimensional manifold where the metric is the push-forward of the logit. There the Riemannian manifold is (0,1) equipped with the metric g<sub>p</sub>(u,v) = uG(p)v where

<footnote n="17">¹⁷A. Hyvärinen and E. Oja, “Independent component analysis: Algorithms and applications,” en, *Neural Networks*, vol. 13, no. 4-5, pp. 411–430, Jun. 2000.</footnote>

<footnote n="18">¹⁸J.-B. Schiratti, S. Allassonnière, O. Colliot, et al., “A Bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations,” *The Journal of Machine Learning Research*, vol. 18, no. 1, pp. 4840–4872, Jan. 2017.</footnote>

<footnote n="19">¹⁹I. Koval, “Learning Multimodal Digital Models of Disease Progression from Longitudinal Data: Methods & Algorithms for the Description, Prediction and Simulation of Alzheimer’s Disease Progression,” en, Ph.D. dissertation, Institut Polytechnique de Paris, Jan. 2020.</footnote>

<footnote n="20">²⁰C. R. Jack, D. S. Knopman, W. J. Jagust, et al., “Update on hypothetical model of Alzheimer’s disease biomarkers,” *Lancet neurology*, vol. 12, no. 2, pp. 207–216, Feb. 2013.</footnote>
$G(p) = \frac{1}{p^2(1-p)^2}$. For a $d$-dimensional dataset, we use the product manifold of this 1-dimensional manifold.

The specific formulation of the curve of one dimension $k$ for one individual $i$ at time $t_{ij}$ in the model is the following:

$$ \eta_k(\psi_i(t_{ij}); \mathbf{w}_i) = \left( 1 + \left( \frac{1}{p_k} - 1 \right) \exp \left( - \frac{v_k(e^{\xi_i}(t_{ij} - t_0 - \tau_i) + t_0) + w_{ik}}{p_k(1 - p_k)} \right) \right)^{-1} \quad (1.1.4) $$

where $(p_k, v_k, t_0)$ are the fixed-effects parameters for the average trajectory on item $k$ defined by $\gamma_k(t_0) = p_k$ and $\dot{\gamma}_k(t_0) = v_k$. $(\tau_i, \xi_i, w_{ik})$ constitute the random effects, $\xi_i$ being the log-acceleration and $\tau_i$ the time-shift both modulating the individual disease timeline while $w_{ik}$ is the space-shift parameter described in the previous section.

In Figure 1.3 we can observe the typical trajectories on the logistic curves manifold in 2 dimensions, as well as the trajectories as a function of time.

### 1.1.3 Bayesian model

The DCM model is a generative model. Therefore there is a probabilistic model on top of the geometrical model. One thing to note is that the DCM model is not structurally identifiable as there can be several triplets $(t_0, \mathbf{p}, \mathbf{v})$ describing the same set of trajectories. Identifiability comes with the probability distributions over the parameters. The authors in$^{21}$ chose the DCM model to be Bayesian, in order to have posterior distributions over the parameters rather than point estimates. We provide next the list of priors. Some new parameters are introduced to have priors over individual parameters.

#### Priors

We first need to introduce new notations to have priors respecting the definition domains of the variables: $p_k = \frac{1}{1+g_k}$, $g_k = e^{\tilde{g}_k}$ and $v_k = e^{\tilde{v}_k}$. This will result in a log-normal prior for $g_k$ and $v_k$, while the other latent parameters will have a normal prior. Here is the complete prior list:

[Population parameters (fixed effects)]

$$ \tilde{g}_k \sim \mathcal{N}(\bar{g}_k, \sigma_g^2) \quad (1.1.5) $$

$$ \tilde{v}_k \sim \mathcal{N}(\bar{v}_k, \sigma_v^2) \quad (1.1.6) $$

$$ \beta_{ml} \sim \mathcal{N}(\bar{\beta}_{ml}, \sigma_{\beta}^2) \quad (1.1.7) $$

[Individual parameters (random effects)]

$$ \tau_i \sim \mathcal{N}(0, \sigma_{\tau}^2) \quad (1.1.8) $$

$$ \xi_i \sim \mathcal{N}(0, \sigma_{\xi}^2) \quad (1.1.9) $$

$$ s_{il} \sim \mathcal{N}(0, 1) \quad (1.1.10) $$

<footnote>$^{21}$J.-B. Schiratti, S. Allassonnière, A. Routier, et al., “A Mixed-Effects Model with Time Reparametrization for Longitudinal Univariate Manifold-Valued Data,” en, in *Information Processing in Medical Imaging*, S. Ourselin, D. C. Alexander, C.-F. Westin, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2015, pp. 564–575.</footnote>
[Noise model]

$$y_{ijk} \sim \mathcal{N}(\eta_k(\psi_i(t_{ij}); \mathbf{w}_i), \sigma_k) \qquad (1.1.11)$$

with many new parameters added. We provide below a summary of those parameters.

## Parameter classes

We divide the parameters into three classes, which is determined by their role in the estimation algorithm, namely the MCMC-SAEM. More details about the algorithm will be provided in the dedicated chapter 5. For now it suffices to understand that it is an iterative algorithm alternating between two steps. In the first step, the expectation (E-step) is approximated stochastically by sampling some variables which are called *latent*. The second step is the maximization (M-step), which consists in maximizing the likelihood over all the non-latent variables, which are called *model parameters*. In the case of the DCM, all the parameters can thus be regrouped in three categories:

* hyperparameters $(\sigma_g, \sigma_v, \sigma_\beta)$ are fixed by the user ahead of estimation. In practical implementation they are all fixed to 0.01 by default. Their role in the prior distribution of $(\tilde{g}, \tilde{v}, \tilde{\beta})$ is to regulate the exploration of the $(\mathbf{g}, \mathbf{v}, \beta)$ parameters during the sampling

* model parameters $\theta = (\tilde{g}, \tilde{v}, \tilde{\beta}, t_0, \sigma_\tau, \sigma_\xi, \sigma_k)$ are all maximized during the M-step in a closed form. Intuitively $(\tilde{g}, \tilde{v}, \tilde{\beta})$ are the mean value of $(\tilde{g}, \tilde{v}, \tilde{\beta})$, and the main purpose of these priors is to control exploration. For individual parameters $(\tau, \xi, \mathbf{s})$ we have a centered gaussian diagonal structure of the random effects. Finally $\sigma_k$ is the standard deviation of the noise on item $k$.

* latent parameters $z = (z_{pop}, (z_i)_{1 \le i \le N})$ are sampled during the E-step. The standard mixed-effect distinction is made between population parameters $z_{pop} = (\tilde{g}, \tilde{v}, \beta)$ and individual parameters $z_i = (\tau_i, \xi_i, \mathbf{s}_i)$.

The implementation of the DCM model can be found in the open-source library Leaspy:

https://leaspy.readthedocs.io/en/stable/

The DCM is still expanding with many new variants, including all the work presented in this manuscript. The Git link to the code is:

https://gitlab.com/icm-institute/aramislab/leaspy.

In terms of vocabulary, we use the term *calibration* to refer to the estimation of the model parameters and latent parameters on a particular dataset. We use the term *personalization* to designate the estimation of only individual parameters with a fixed set of model parameters and population parameters. Typically, we use a training dataset to calibrate a DCM model, and then we use a test set for predictions. The test set is composed only of individuals that were not used during calibration, therefore we have no individual parameters for them, hence the need for the personalization. We fix the parameters of the DCM that was already calibrated, and we estimate the best individual parameters for the patients in the test set. Finally we can perform predictions, imputation, clustering or covariate analysis for these new patients.
1.1.4 Limitations

One could argue that the logistic curves model is sufficient for the needs of most neurodegenerative applications as showcased by the results obtained in Alzheimer's disease (AD)<fnref n="22" />. In this work the logistic curves DCM outperformed 56 other prediction methods in the TADPOLE challenge, where the goal was to predict the changes of cognitive scores and biomarkers for AD patients several years in advance. Other successful applications have been made in Huntington's disease<fnref n="23" /> and Parkinson's disease<fnref n="24" />.

However this approach is still very rigid in the sense that we hypothesise the shape of the trajectories. The logistic shape can be a byproduct of a biological process model, which is the case for the accumulation of prions in the brain for instance, where an accurate mathematical model of prion spread in the cortex<fnref n="25" /> shows that the overall quantity of misfolded proteins evolves closely to a logistic over time. For cognitive scores, the logistic also unfolds from the design of the score. Indeed they are often created such that patients are normally distributed across its scale, leading to the cumulative distribution function being akin to a logistic. But for other markers the logistic hypothesis is not properly grounded, as is the case for imaging and cerebrospinal fluid markers in AD. The model would thus benefit from a wider range of possible trajectories.

Another limitation comes from the injectivity and the product manifold. Since the model is injective and is the result of a *d*-dimensional product of the 1-dimensional setting, the trajectory has to be injective in each dimension. This results in the trajectories being monotonous, due to their continuity. This means that we can only select markers with a monotonous progression over time. Even though usually this is no problem with a slow progressive incurable disease, some markers that are relevant to the disease might have fluctuating levels.

1.2 Geodesic Bending

These limitations led to a work questioning the flexibility of the model, and whether we could try
to learn the metric instead of imposing it. The proposed solution for this metric learning problem
led to the model proposed in this section. The method was developed in collaboration with Samuel
Gruffaz, which led to the paper<fnref n="26" /> from which comes the materials of the current chapter. My main
contribution to this work is in the estimation algorithm, which is described in section 5.2.2. The
conception of the model and the experiments as well as the analysis of the results was an equal
contribution of both of us, while the practical experiments were carried out by Samuel Gruffaz.
The mathematical properties of the model can be hard to understand without some notions of
Riemannian geometry, so we encourage the reader to have a look at the dedicated appendix section
15.

<footnote n="22">²²I. Koval, A. Bône, M. Louis, et al., "AD Course Map charts Alzheimer's disease progression," eng, *Scientific Reports*, vol. 11, no. 1, p. 8020, Apr. 2021.</footnote>

<footnote n="23">²³I. Koval, T. Dighiero-Brecht, A. J. Tobin, et al., "Forecasting individual progression trajectories in Huntington disease enables more powered clinical trials," en, *Scientific Reports*, vol. 12, no. 1, p. 18928, Nov. 2022.</footnote>

<footnote n="24">²⁴R. Couronné, A. Valladier, M. Vidhaillet, et al., "Modeling the progression of Parkinson's Disease: Comparison of subjects with and without Sleep Disorders," en, p. 4,</footnote>

<footnote n="25">²⁵J. Weickenmeier, M. Jucker, A. Goriely, et al., "A physics-based model explains the prion-like features of neu-
rodegeneration in Alzheimer's disease, Parkinson's disease, and amyotrophic lateral sclerosis," en, *Journal of the*
*Mechanics and Physics of Solids*, vol. 124, pp. 264–281, Mar. 2019.</footnote>

<footnote n="26">²⁶S. Gruffaz, P.-E. Poulet, E. Maheux, et al., “Learning Riemannian metric for disease progression modeling,” in *Advances in Neural Information Processing Systems*, M. Ranzato, A. Beygelzimer, Y. Dauphin, et al., Eds., vol. 34, Curran Associates, Inc., 2021, pp. 23780–23792.</footnote>
### 1.2.1 Introduction and related work

As debated previously, the rigid setting of the DCM model comes from the metric choice which induces the shape of the trajectory curves. The two main use cases of the DCM model are the linear case and the logistic curves model. But what if we could learn the metric instead of fixing it? This has been attempted before with the DCM model. The first approach for metric learning was proposed in<fnref n="27" /> by learning of the metric directly as parameters of the model and later extended in<fnref n="28" />. The application of the DCM to MRI images, where the manifold cannot be specified, has led to the use of a variational auto-encoder (VAE) to learn the metric implicitly in<fnref n="29" />. In this last paper, individual trajectories were seen as transformations of straight lines evolving in a latent space learned by the auto-encoder. Although this semi-parametric representation enabled multiple progression profiles, it was not designed for disentangling space and time variability. On top of that, both approaches were computationally heavy and estimation was burdened by a high number of parameters requiring sampling.

A standard method for trajectory learning is to use the theory of shape analysis and more especially the deformation of shapes in time. Authors commonly learn the deformation as a diffeomorphism with a Reproducing Kernel Hilbert Space (RKHS), the kernel regularity encoding the smoothness of the deformation through time<fnref n="30" />,<fnref n="31" />. We believe that this deformation method is an interesting alternative to the auto-encoders to increase model flexibility while keeping control on space and time variability. Another approach would be to rely on a push-forward method to learn a metric with pseudo-diffeomorphic transformations<fnref n="32" />.

In this work, following<fnref n="33" /> but taking a step back, we propose a semi-parametric method using a RKHS to learn the Riemannian metric of DCM models. We thereby retain the possibility to learn the inter-variability of patient trajectories thanks to the mixed-effect framework, keeping the advantages of the geometric and Bayesian approach presented in<fnref n="34" />. First, we present the model for learning the metric. The modifications of the estimation algorithm of the DCM are tackled in the dedicated chapter 5.2.2. We validate the presented method on synthetic data and on a real dataset by comparing it with previous models on the task of predicting Alzheimer's disease patient's biomarker progression. Finally, we discuss limitations and possible future works in the discussion section.

<footnote n="27">²⁷M. Louis, “Computational and statistical methods for trajectory analysis in a Riemannian geometry setting,” en, Ph.D. dissertation, Sorbonne Université, Oct. 2019.</footnote>

<footnote n="28">²⁸B. Sauty and S. Durrleman, “Riemannian Metric Learning for Progression Modeling of Longitudinal Datasets,” in *2022 IEEE 19th International Symposium on Biomedical Imaging (ISBI)*, Mar. 2022, pp. 1–5.</footnote>

<footnote n="29">²⁹B. Sauty and S. Durrleman, “Progression Models for Imaging Data with Longitudinal Variational Auto-Encoders,” en, in *Medical Image Computing and Computer Assisted Intervention – MICCAI 2022*, L. Wang, Q. Dou, P. T. Fletcher, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer Nature Switzerland, 2022, pp. 3–13.</footnote>

<footnote n="30">³⁰A. Bône, M. Louis, B. Martin, et al., “Deformetrica 4: An Open-Source Software for Statistical Shape Analysis,” Sep. 2018.</footnote>

<footnote n="31">³¹S. Joshi and M. Miller, “Landmark matching via large deformation diffeomorphisms,” *IEEE Transactions on Image Processing*, vol. 9, no. 8, pp. 1357–1370, Aug. 2000.</footnote>

<footnote n="32">³²G. Lebanon, *Learning Riemannian Metrics*, Oct. 2012.</footnote>

<footnote n="33">³³M. Louis, “Computational and statistical methods for trajectory analysis in a Riemannian geometry setting,” en, Ph.D. dissertation, Sorbonne Université, Oct. 2019.</footnote>

<footnote n="34">³⁴J.-B. Schiratti, S. Allassonnière, O. Colliot, et al., “A Bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations,” *The Journal of Machine Learning Research*, vol. 18, no. 1, pp. 4840–4872, Jan. 2017.</footnote>
1.2.2 Method

Notations in all the following are: |·| for the usual Euclidean norm, 〈·, ·〉 the Euclidean scalar product, |||L||| = sup$_{x}$$_{∈}$$_{(}$$_{R}$ $^{d}$) ∗ |L(x)| |x| for all linear function L and |||df |||$_{∞}$ = sup$_{x}$ |||df (x)||| for all differentiable function f .

In this part, M is an open subset of R d equipped with the Riemannian metric g . We highlight objects depending on g with the notation g because g will come to be seen as a variable in the following. For instance we note the reference geodesic γ g $_{0}$, and an individual trajectory γ g $_{i}$.

Traditionally, the metric g is selected using empirical arguments. This metric g is the most important since the whole model depends on it: it encodes the shape of geodesics with g and the variability due to exp-parallelelisation with η $^{g}$. One of the main issues of setting the metric for a 1- dimensional manifold and then taking the product is that product metrics assume that the geodesics coordinates are independent which is not necessarily true when they describe neurodegenerative bio-markers progression$^{35}$.

Thus, we assume that the noise effect ϵ$_{ijk}$ = y$_{ijk}$ − η$_{k}$ ( ψ$_{i}$ ( t$_{ij}$ ); w$_{i}$ ) is in reality different from the real data noise ϵ noise ijk because of the miss-specification of the metric g , in such a way that ϵ$_{ijk}$ = ϵ g ijk + ϵ noise ijk with ϵ g ijk the error due to the modelling. The goal is to minimize ϵ g ijk .

To state the problem, we suppose that the population and individual parameters θ = ( θ, ( z$_{i}$ )$_{i}$) are estimated and we want to minimize the squared reconstruction error L ( g ) according to g :

$$
L(g) := \sum_{ijk} |\epsilon_{ijk}|^2 = \sum_{ijk} |\epsilon_{ijk}^g + \epsilon_{ijk}^{\text{noise}}| = \sum_{ij} |y_{ij} - \eta_k(\psi_i(t_{ij}); \mathbf{w}_i)|^2 \quad (1.2.1)
$$

To prevent the risk of overfitting since ϵ noise ijk is present in L, we can add a regularizing term or constrain the set of admissible metrics g . Next, we describe a method to solve this optimization problem.

The task of learning a Riemannian metric has already been addressed in the literature$^{36,37}$ ,38 and most of the time, it is reduced to the task of learning a diffeomorphism using the following proposition:

Proposition 1.1 (Pushforward metrics) . Provided ( M , g ) a Riemannian space , N a manifold and φ : M → N a C 1 diffeomorphism , we can equip N with the Riemannian metric g φ defined as:

$$
\forall p \in N, \quad \forall w, v \in TN_p, \quad g_p^\phi(w, v) = g_{\phi^{-1}(p)}(d\phi^{-1}(p).w, d\phi^{-1}(p).v)
$$

Moreover, φ is an isometry, which implies that for all ( m, v ) in T M and γ : ( − 1 , 1 ) → M a differentiable curve:

$$
\mathrm{Exp}_{\phi(m)}^{g}(d\phi(m).v) = \phi \circ \mathrm{Exp}_{m}^{g}(v), \quad \mathrm{P}_{\phi \circ \gamma, t_0, s}^{g}(d\phi(\gamma(t_0)).v) = d\phi(\gamma(s)).\mathrm{P}_{\gamma, t_0, s}^{g}(v)
$$

<footnote>$^{35}$B. O. Taddé, H. Jacqmin-Gadda, J.-F. Dartigues, et al. , “Dynamic modeling of multivariate dimensions and their temporal relationships using latent processes: Application to Alzheimer’s disease,” en, Biometrics , vol. 76, no. 3, pp. 886–899, 2020.</footnote>

<footnote>$^{36}$G. Lebanon, Learning Riemannian Metrics , Oct. 2012.</footnote>

<footnote>$^{37}$M. Louis, “Computational and statistical methods for trajectory analysis in a Riemannian geometry setting,” en, Ph.D. dissertation, Sorbonne Université, Oct. 2019.</footnote>

<footnote>$^{38}$S. r. Hauberg, O. Freifeld, and M. Black, “A Geometric take on Metric Learning,” in Advances in Neural Information Processing Systems , vol. 25, Curran Associates, Inc., 2012.</footnote>
Thanks to this property, we can start from an initial metric $g$ in $\mathcal{M}$ and reach a wide variety of Riemannian metrics by considering $(\phi(\mathcal{M}), g_0^\phi)$ with $\phi$ in $\mathcal{F}$ a set of $C^1$-diffeomorphism. Moreover, the Exp-parallelelisation $\eta^{g^\phi}$ can be expressed in closed form as soon as $\phi$ and $\eta^g$ are expressed in closed form as well. Nevertheless, it restricts the possibilities because the curvature is preserved. For example if $g$ is the Euclidean metric and $\mathcal{M} = \mathbb{R}^d$, the space $(\phi(\mathcal{M}), g^\phi)$ will be always flat (curvature equals zero) since $(\mathcal{M}, g)$ is flat.

Learning a diffeomorphism is a more common task than learning a metric especially in the field of shape analysis with the LDDMM algorithm$^{39}$ and even in deep learning with the invertible networks$^{40}$ for different applications$^{41,42}$. Inspired from the LDDMM framework, we chose to use the following proposition to construct the set of $C^1$-diffeomorphism $\mathcal{F}$:

**Proposition 1.2** (Sufficient condition for a diffeomorphism). If $\phi = id + f$, with $f$ a bounded function in $C^1(\mathbb{R}^d, \mathbb{R}^d)$ such that $|||df|||_\infty < 1$ (H1), then $\phi$ is a $C^1$ diffeomorphism.

In all the following we choose $\mathcal{M} = \mathbb{R}^d$ which is the most flexible choice. For example if we want to subsequently work on $(0,1)$, consider $\phi' = \phi \circ \sigma$ with $\sigma = \frac{1}{1+\exp(-t)}$. To respect the condition of the previous proposition, we choose to take $f$ in an RKHS denoted $\mathcal{H}$ to encode its regularity in its norm. Moreover, the RKHS functions are handy in non-parametric optimization thank to the representer theorem (for a comprehensive review of kernel methods$^{43}$). We derived the following conditions on the kernels to fit the constraints:

**Lemma 1.1.** If $k$ is a kernel such that $k(x, y) = g(x-y)$ Id with $g \in C^2(\mathbb{R}^d, \mathbb{R})$ and $h(x, y) = g(x-y)$ is a bounded kernel, then $\forall f \in \mathcal{H}$ s.t. $||f||_\mathcal{H} < \frac{1}{\sqrt{|\nabla g(0)|}}$, $|||df|||_\infty < 1$ and $f$ is bounded (H1).
More specifically for the common kernels:

* The Gaussien kernel; $g(x) := \exp(-\frac{|x|^2}{2\sigma^2})$, $\sigma > 0$, $\frac{1}{|\nabla g(0)|} = \frac{\sigma}{\sqrt{d}}$.

* The Sobolev kernel or generalized T-student kernel, $g(x) = \frac{1}{(1+\frac{|x|}{2\sigma^2})^a}$, $\sigma > 0, a > d$, $\frac{1}{|\nabla g(0)|} = \frac{\sigma}{a\sqrt{d}}$

Proof.

$$ \begin{align*} \partial_i f(x).\alpha &= \partial_i k(x, \cdot)\alpha f \\ df(x)^T.\alpha &= (\partial_i k(x, \cdot)\alpha f)_i \\ ||df(x)^T.\alpha||_2^2 &= \sum_{i=1}^d \partial_i k(x, \cdot)\alpha f^2 \\ \forall \alpha \in S^{d-1}, \quad ||df(x)^T.\alpha||_2^2 &\le ||f||_H^2 \sum_{i=1}^d ||\partial_i k(x, \cdot)\alpha||_H^2 \end{align*} $$

<footnote>$^{39}$S. Joshi and M. Miller, “Landmark matching via large deformation diffeomorphisms,” *IEEE Transactions on Image Processing*, vol. 9, no. 8, pp. 1357–1370, Aug. 2000.</footnote>

<footnote>$^{40}$J.-H. Jacobsen, A. Smeulders, and E. Oyallon, *I-RevNet: Deep Invertible Networks*, Feb. 2018.</footnote>

<footnote>$^{41}$L. Younes, “Diffeomorphic Learning,” *Journal of Machine Learning Research*, vol. 21, no. 220, pp. 1–28, 2020.</footnote>

<footnote>$^{42}$M. A. Rana, A. Li, D. Fox, et al., “Euclideanizing Flows: Diffeomorphic Reduction for Learning Stable Dynamical Systems,” en, in *Proceedings of the 2nd Conference on Learning for Dynamics and Control*, PMLR, Jul. 2020, pp. 630–639.</footnote>

<footnote>$^{43}$T. Hofmann, B. Schölkopf, and A. J. Smola, “Kernel methods in machine learning,” *The Annals of Statistics*, vol. 36, no. 3, pp. 1171–1220, Jun. 2008.</footnote>
$$ \forall \alpha \in S^{d-1}, |||df(x)|||^2 \le ||f||_H^2 \sum_{i=1}^{d} ||\partial_i k(x, \cdot)\alpha||_H^2 $$

We can enforce $||f||_H^2 \sum_{i=1}^d ||\partial_i k(x, \cdot)\alpha||_H^2 < 1$ by assuming:

$$ \forall \alpha \in S^{d-1}, ||f||_H^2 < \frac{1}{\sum_{i=1}^d ||\partial_i k(x, \cdot)\alpha||_H^2} $$

Using the properties of a $C^2$ kernel, we have:

$$ ||\partial_i k(x, \cdot)\alpha||_H^2 = \alpha^T \partial_{i,i}^2 k(x, \cdot)\alpha $$

If $k(x, y) = h(x, y)I_d$, we have:

$$ \alpha^T \partial_i^2 k(x, x)\alpha = ||\alpha||^2 \partial_{i,i}^2 h(x, x) = \partial_{i,i}^2 h(x, x) $$

Then

$$ ||f||_H^2 < \frac{1}{\Delta h(x, x)} $$

With $h(x, y) = g(x - y)$, we have

$$ ||f||_H^2 < \frac{1}{\Delta g(0)} $$

□

In all the following we choose a kernel $k$ respecting the previous conditions such that for all $f$ in the closed ball $\bar{B}_H(0, c)$ with $c$ a positive constant, (H1) is verified. Note that we can prevent overfitting thanks to this constraint.

## Estimation overview

If the population and random effects $(z_{pop}, (z_i))$ are estimated, the minimization of the error of reconstruction over the metric $g$ can be performed by solving the following problem:

$$ f^* \in \underset{f \in \bar{B}_H(0, c)}{\operatorname{argmin}} L(f), L(f) = \sum_{i,j} |y_{i,j} - \gamma_i(t_{i,j}) - f(\gamma_i(t_{i,j}))|^2 \quad (1.2.2) $$

With the representer theorem of the RKHS, we have:

$$ \forall x \in f^*(x) := \sum_{i,j} k(x_{i,j}, x) w_{i,j}, x_{i,j} = \gamma_i(t_{i,j}) \quad (1.2.3) $$

This allows us to reduce the metric optimization problem to a convex optimization problem, which can be solved. The number of weights in the sum might be too large and cause overfitting, so we restrain the basis functions to a lower number of $N_c$ centered on control points. More details about this in the section about optimization 5.2.2. In the end we obtain a new algorithm called alternating maximization. It alternates between several iterations of the MCMC-SAEM algorithm to optimize the DCM parameters and a metric optimization. This algorithm can be performed several times in a row to compose diffeomorphisms. We start with a simple metric $g_{init}$ (usually the Euclidean metric or the metric resulting from the push-forward of the logit), and everytime we
apply the algorithm amounts to a composition of the previously estimated diffeomorphism $\phi$ with the newly estimated diffeomorphism.

For $N_{\text{comp}}$ compositions, we have $g = g^{\Phi_{N_{\text{comp}}}}$ with $\Phi_{N_{\text{comp}}} = \phi_{N_{\text{comp}}} \circ \dots \circ \phi_1$ where $\phi_i = id + \sum_{j=1}^{N_c} k(x_j^i, \dots, w_j^i, (x_j^i)_j)$ the $N_c$ control points and $(w_j^i)_j$ their associated weights. This structure resembles the deep neural networks considering $N_c$ as the width and $N_{\text{comp}}$ as the depth. The total time complexity is $O(n_{\text{MCMC}}dN_cN_{\text{comp}}^3)$ which reduces in practice the choice of $N_{\text{comp}}$. In theory, the runtime is $N_{\text{comp}}$ times the normal runtime of a standard DCM, but this can be improved in practice. We called this method Geodesics Bending (GB).

### 1.2.3 Experiments

All the methods are developed in Python by extending the open-source Leaspy library: https://leaspy.readthedocs.io and run on a 2.80GHz CPU with 16 GB RAM.

#### Synthetic data

In this part, we study the learning capacity of the method and its stability on synthetic data depending on the values of the depth $N_{\text{comp}}$, the width $N_c$ and the initial metric $g_{\text{init}}$. We study the algorithm in the context of neurodegenerative diseases. The parameters are selected to be realistic.

In all the following, to perform the metric estimation, we choose $k$ to be the Gaussian kernel (experiments with the Sobolev kernel gives similar results). The Gaussian kernel is parametrized by its bandwidth $\sigma$, which is an hyperparameter that requires tuning. $\sigma$ controls the size of the radial basis functions, so a small $\sigma$ allows for sharp local variations in the diffeomorphism while a large $\sigma$ only allows slow variations. Then we chose $n_{\text{MCMC}} = 200$ MCMC-SAEM iterations between metric optimization steps ($n_{\text{MCMC}} = 10000$ at the very first step since the mixed-effects are not necessarily well initialized whereas for subsequent iterations the DCM parameters starts closer to the optimum), $\theta_{\text{init}}$ is computed empirically from patients features$^{44}$. It leaves $N_{\text{comp}}$, $N_c$ and $g_{\text{init}}$ to select according to the situation.

**Data generation** Once we have selected the number of patients $N_{\text{pat}}$ and the observations' dimension $d$, we generate data according to the DCM model selecting the fixed effects:

$$t_0, \mathbf{v}, \mathbf{p}, g_{\text{gen}}, \sigma_{\tau}^2, \sigma_{\xi}^2, \sigma_{\text{noise}}^2$$

and sampling the random effects parameters such that:

$$\tau_i \sim \mathcal{N}(0, \sigma_{\tau}^2), \xi_i \sim \mathcal{N}(0, \sigma_{\xi}^2), \epsilon_{i,j} \sim \mathcal{N}(0, \sigma_{\text{noise}}^2), s_i^m \sim \mathcal{N}(0, 1), w_i = \sum_m s_i^m A_m$$

and time parameters $\sigma_t$, $n_{\text{min}}^t$, $n_{\text{delta}}^t$ to generate the patient time-points $(t_{i,j})_{1 \le j \le T_i}$ randomly:

$$T_i = n_{\text{min}}^t + \delta_i, \delta_i \sim B(n_{\delta}, \frac{1}{2}), t_{i,j} \sim \mathcal{N}(t_0, \sigma_t), y_{i,j} = \gamma_i^{g_{\phi}}(t_{i,j}) + \epsilon_{i,j}$$

For both experiments, we set $N_{\text{pat}} = 500$, $d = 2$, $n_{\text{min}}^t = 2$, $n_{\text{delta}}^t = 6$, $\sigma_t = 4$.

<footnote>$^{44}$I. Koval, “Learning Multimodal Digital Models of Disease Progression from Longitudinal Data: Methods & Algorithms for the Description, Prediction and Simulation of Alzheimer’s Disease Progression,” en, Ph.D. dissertation, Institut Polytechnique de Paris, Jan. 2020.</footnote>
**Convergence** First, we experiment whether learning the metrics with kernels enables the regression of exotic trajectories. We choose $\sigma_{\text{noise}} = 0.01$, $g_{\text{gen}} = g^{\phi}$ with $\phi(x) := \exp(0.07(x+\sin(x))-1)$ and $g = \langle \cdot, \cdot \rangle$ for the data generation and $N_c(\sigma = 0.08) \approx 90$ for the metric estimation.

<figure><img src="image_13.png" /><figcaption>Figure 1.4: (a) The average trajectory on the second coordinate is displayed for each value of $N_{\text{comp}}$ in the Convergence experiment. (b) The average trajectory on the second coordinate is displayed for each value of $N_{\text{comp}}$ in the Stability experiment. Each step corresponds to a new diffeomorphism composition.</figcaption></figure>

Increasing the number of compositions improves the capacity to recover the true average trajectory as we can see in Figure 1.4. The estimated noise variance $\tilde{\sigma}_{\text{noise}}^2$ decreases linearly and plateaus around the true value of noise $\sigma_{\text{noise}}^2 = 0.013$ which shows that the model is a little under-fitting. The kernel regularization impacts the smoothness of the function, making it harder to fit the fast oscillations. By testing different values of the width $N_c$, we observed that the lower the kernel variance $\sigma$, the higher the risk to overfit. Conversely, the higher $\sigma$, the higher the risk to underfit. As for deep neural networks, there is an optimum to find regarding the hyper-parameters. In practise with datasets on neurodegenerative disease, normalized data are very noisy $\sigma_{\text{noise}}^2 \approx 0.05$. To avoid over-fitting, we select the kernel variance $\sigma$ in [0.1,0.5]. The choice of $g_{\text{init}}$ influences the number of compositions required to reach a good approximation of the "true" metric.

**Stability** Secondly, we experiment whether the method is stable: beginning from the metric which has generated the data with $g_{\text{init}}$, we observe the effect of supernumerary compositions. We choose $g_{\text{gen}} = g^{\text{sig}}$ with $\text{sig}(x) := \frac{1}{1+\exp(-x)}$ and $\sigma_{\text{noise}} = 0.05$ for the data generation, $g_{\text{init}} = g^{\text{sig}}$ and $N_c(\sigma = 0.2) \approx 31$ for the metric estimation.

The compositions on $\phi$ produce local fluctuation on the sigmoid trend as pictured in figure 1.4. There are less observations near the inflexion point which causes over-fitting. This phenomenon is also present in datasets on neurodegenerative disease: areas of the trajectory where data are scarce can be over-fitted by GB. Nevertheless, the value of the DCM parameters are nearly constant after its first estimation in the algorithm which strengthens the model interpretability and encourages the reduction of the number of MCMC-SAEM steps $n_{\text{MCMC}}$.
**The problem of generalization** The method reveals to be flexible and quite stable provided the hyper-parameters are wisely selected. To assess whether the learning complexity of the method makes sense in practice, we propose to measure its capacity of generalization by recording its performance on prediction tasks within a 5-folds cross validation framework: running the algorithm on 4-folds, predicting on the last fold the last patient visit from its first visits, or predicting the previous visits from the last visits (as for data imputation). In the next part, these experiences are carried out on a reference dataset of neurodegenerative diseases with different types of data (cognitive scores, sub-cortical volumes from MRI, cerebro-spinal fluid biomarkers) to show the potential of a flexible method such as the one we introduced.

## Real data

**Introduction to ADNI** The data used in this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu), launched in 2003 as a public-private partnership, led by Principal Investigator Michael W. Weiner, MD. The primary goal of ADNI has been to test whether serial magnetic resonance imaging (MRI), positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimer's disease (AD).

First, the algorithm was tested on the subset of AD patients. For each visit, we jointly model their Mini-Mental State Examination (MMSE, a cognitive score), the volume of their hippocampi and lateral ventricles normalized by their intracranial volume (HIPP and VENTS respectively) and the concentration of amyloid-β$_{1-42}$ in their cerebro-spinal fluid (ABETA). The goal was to understand the impact of the metric estimation on these typical AD biomarkers. Then, the model was compared to the state-of-the-art with The Alzheimer's Disease Prediction Of Longitudinal Evolution (TADPOLE) challenge$^{45}$.

**An Alzheimer's disease cohort** The cohort is composed of patients from all ADNI phases (ADNI1, ADNI GO, ADNI2, ADNI3), having been diagnosed with AD with a minimum of 2 visits (4.15 visits on average and 1.5 visits for standard deviation) and $ABETA < 977pg/mL$ cutoff (so to get amyloid-positive patients only$^{46,47}$). The dataset statistics for the selected variables are shown in Table 1.1.

Sub-cortical volumetric segmentations of MRI were performed with the Freesurfer image analysis suite v6.0.0, which is documented$^{48}$ and freely available for download online: http://surfer.nmr.mgh.harvard.edu/. All scores were normalized in [0,1] (to motivate the logistic prior for trajectories), using preprocess-ing tools from the scikit-learn library$^{49}$: MMSE was affinely-mapped using the score bounds, other

<footnote>$^{45}$R. V. Marinescu, N. P. Oxtoby, A. L. Young, et al., "TADPOLE Challenge: Accurate Alzheimer's disease prediction through crowdsourced forecasting of future data," eng, *PRedictive Intelligence in MEdicine. PRIME (Workshop)*, vol. 11843, pp. 1-10, Oct. 2019.</footnote>

<footnote>$^{46}$C. R. Jack, D. A. Bennett, K. Blennow, et al., "A/T/N: An unbiased descriptive classification scheme for Alzheimer disease biomarkers," eng, *Neurology*, vol. 87, no. 5, pp. 539-547, Aug. 2016.</footnote>

<footnote>$^{47}$O. Hansson, J. Seibyl, E. Stomrud, et al., "CSF biomarkers of Alzheimer's disease concord with amyloid-PET and predict clinical progression: A study of fully automated immunoassays in BioFINDER and ADNI cohorts," eng, *Alzheimer's & Dementia: The Journal of the Alzheimer's Association*, vol. 14, no. 11, pp. 1470-1481, Nov. 2018.</footnote>

<footnote>$^{48}$B. Fischl, D. H. Salat, E. Busa, et al., "Whole brain segmentation: Automated labeling of neuroanatomical structures in the human brain," eng, *Neuron*, vol. 33, no. 3, pp. 341-355, Jan. 2002.</footnote>

<footnote>$^{49}$F. Pedregosa, G. Varoquaux, A. Gramfort, et al., "Scikit-learn: Machine Learning in Python," en, *MACHINE LEARNING IN PYTHON*,</footnote>
<table><caption>Table 1.1: Statistics for ADNI experiment. MMSE: Mini-mental state examination. MRI X ICV: intracranial volume of region X observed on MRI.</caption><thead><tr><th></th><th>MMSE</th><th>β-amyloïd 42</th><th>MRI Hippocampus ICV</th><th>MRI Ventricles ICV</th></tr></thead><tbody><tr><td>Count</td><td>3799</td><td>1019</td><td>3395</td><td>3392</td></tr><tr><td>Mean</td><td>23.78</td><td>703.6</td><td>0.00379</td><td>0.0321</td></tr><tr><td>Standard deviation</td><td>4.70</td><td>385.6</td><td>0.00072</td><td>0.0141</td></tr><tr><td>Min</td><td>0</td><td>0</td><td>0.001</td><td>0.0053</td></tr><tr><td>25%</td><td>22</td><td>477.1</td><td>0.00327</td><td>0.0217</td></tr><tr><td>50%</td><td>25</td><td>612.4</td><td>0.00373</td><td>0.0301</td></tr><tr><td>75%</td><td>27</td><td>782.0</td><td>0.00425</td><td>0.0401</td></tr><tr><td>Max</td><td>30</td><td>3139.0</td><td>0.00726</td><td>0.0903</td></tr></tbody></table>

biomarkers were Box-Cox transformed<fnref n="50" /> to un-skew their distribution, clipped between their 0.1 and 99.9 percentiles and finally affinely-mapped into [0,1]. When needed, scores were reversed such that they all increase during disease progression. The models being generative, we did not drop out missing values. We compare the DCM models using the Riemannian metric $g^{\text{sig}}$ (baseline) and GB ($g_{\text{init}} = g^{\text{sig}}$, $N_{\text{comp}} = 6$ and $N_c(\sigma = 0.24)$), both using 2 principal directions of variation for space-shifts. Note that the baseline (DCM in yellow in Figure 1.6) is the very first iteration of GB optimization, before any diffeomorphism is applied.

Looking at the Figure 1.6, the MMSE trajectory learnt by GB has a more exponential shape compared to the DCM and the hippocampus volume trajectory has been transformed from a nearly linear curve into a piece-wise linear curve, the inter-patient variability is changed correspondingly (see Figure 1.5). It shows that there are different trajectories behaviors that DCM cannot take into account. With GB, at the beginning of the disease we observe that hippocampus volume decreases while MMSE stays constant, which is clinically coherent. We observe that the algorithm begins to stabilize after 4 steps which motivates that $N_{\text{comp}} = 6$ is a good trade-off between under-fitting and over-fitting.

**Experiments: (1) Future visits prediction and (2) data imputation.** For both experiments, models are trained with a 5-folds cross-validation. In the prediction experiment, we estimate, in the out-of-train fold, the individual parameters $z_i$ on the first $n_f$ visits of patients and predict on their very last visit. In the imputation experiment, in the out-of-train fold, we estimate the individual parameters $z_i$ on all the patients visits except $n_r$ visits drawn uniformly in its set of observations and predict on these $n_r$ removed visits. We evaluate the performance of estimations with the Mean Absolute Errors (MAE) on the 5 test-sets, depending on the value of $n_r$ and $n_f$ varying from 1 to 2.

**Results.** The table 1.2 demonstrates that GB outperforms DCM except on $\beta$-amyloïd, no matter the value of $n_f$ and $n_r$.

However, GB fails to generalize on $\beta$-amyloïd because of the small number of data (576 missing values out of 909 visits). When few information is available logistic prior for trajectories' shape is already almost optimal. It has been remarked that GB makes better predictions than DCM especially on patients with high scores, this fact is worth mentioning if we want to later improve

<footnote n="50">⁵⁰G. E. P. Box and D. R. Cox, “An Analysis of Transformations,” en,</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 1.5: Estimated tubular coordinate system for the hippocampus volume and MMSE scores. This shows the distribution of the spatial variability around the average geodesics in red by changing the space-shifts on a the first principal direction ($w_i = sA_1, s \in [-1, 1]$) in plots **a** and **b**, and on the second direction ($w_i = sA_2, s \in [-1, 1]$) in plots **c** and **d**. 60 and 87 years are the 0.05 and 0.95 quantile for the ages ($t_{i,j}$). Notice how the second source allowed the GB model to be non-injective on the MMSE (plot **d**).</figcaption></figure>

predictions with ensemble methods<fnref n="51" />. Now, we should assess the proposed method on a more general cohort including MCI and cognitively normal patients.

**TADPOLE Challenge.**<fnref n="52" /> The TADPOLE training set is composed of data from the first three ADNI phases (ADNI 1, ADNI GO and ADNI 2). It includes approximately 1500 features acquired

<footnote n="51">⁵¹T. G. Dietterich, “Ensemble Methods in Machine Learning,” en, in *Multiple Classifier Systems*, ser. Lecture Notes in Computer Science, Berlin, Heidelberg: Springer, 2000, pp. 1–15.</footnote>

<footnote n="52">⁵²R. V. Marinescu, N. P. Oxtoby, A. L. Young, et al., “TADPOLE Challenge: Accurate Alzheimer’s disease prediction through crowdsourced forecasting of future data,” eng, *PRedictive Intelligence in MEdicine. PRIME (Workshop)*, vol. 11843, pp. 1–10, Oct. 2019.</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 1.6: Compositions effect on the learning of each marker. The effects of each learning step on the average trajectory $\gamma_0$ is shown.</figcaption></figure>

<table><caption>Table 1.2: Experiment performances recorded with the MAE. The proposed method is compared with the alternative using a paired, two-sided Wilcoxon signed rank test. *: $p < 0.05$, **: $p < 0.01$, ***: $p < 0.001$.</caption><thead><tr><th rowspan="2">Experiment</th><th colspan="2">MMSE</th><th colspan="2">HIPP</th><th colspan="2">VENTS</th><th colspan="2">ABETA</th></tr><tr><th>GB</th><th>DCM</th><th>GB</th><th>DCM</th><th>GB</th><th>DCM</th><th>GB</th><th>DCM</th></tr></thead><tbody><tr><td>Pred(n<sub>f</sub> = 1):</td><td>0.140</td><td>0.146</td><td><strong>0.046</strong>***</td><td>0.056</td><td><strong>0.042</strong>***</td><td>0.056</td><td>0.124</td><td>0.114</td></tr><tr><td>Pred(n<sub>f</sub> = 2):</td><td><strong>0.132</strong>**</td><td>0.146</td><td><strong>0.037</strong>***</td><td>0.046</td><td><strong>0.038</strong>***</td><td>0.049</td><td>0.123</td><td><strong>0.109</strong>*</td></tr><tr><td>Imp(n<sub>r</sub> = 1):</td><td><strong>0.080</strong>**</td><td>0.088</td><td><strong>0.029</strong>***</td><td>0.31</td><td><strong>0.017</strong>***</td><td>0.020</td><td>0.128</td><td>0.123</td></tr><tr><td>Imp(n<sub>r</sub> = 2):</td><td><strong>0.081</strong>***</td><td>0.089</td><td><strong>0.028</strong>***</td><td>0.030</td><td><strong>0.023</strong>***</td><td>0.027</td><td>0.120</td><td><strong>0.111</strong>**</td></tr></tbody></table>
from 1737 subjects (957 males and 780 females) during 12.741 visits for at most 22 distinct time points, between 2003 and 2017. Challengers were evaluated on their prospective predictions on enrolled individuals rolling-over for ADNI next phase, regarding the Alzheimer's Disease Assessment Scale Cognitive Score (ADAS-Cog13) and the volume of lateral ventricles normalized by the intracranial volume (VENTS), with the MAE metric. The final test set, disclosed in 2019 (time to prediction is about 2 years), is composed of 223 follow-up visits. This challenge allows for a fair comparison of the prediction performances presented here against 56 alternative methods. These methods include penalized regression, linear mixed-effect models, recurrent neural networks, and multi-task learning. Since the challenge was closed, more algorithms have been presented for the same prediction task, using the same dataset, including<sup>53</sup> and<sup>54</sup> using deep RNNs, and<sup>55</sup> using random forests. The forecast of ventricle volumes (VENTS) and cognitive decline (ADAS-cog) presented in this work have smaller mean absolute errors than all these competing methods.

We compared GB and DCM to the best Challenger. For the training, we select all patients having at least 2 visits and being amyloid-positive (ABETA < 977pg/mL<sup>56</sup>), it left us 765 patients with 5.5 visits on average. We use the observations of ADAS-Cog13, MMSE, hippocampus volume, ventricles volume and CSF P-TAU to calibrate the models with a 5-folds cross-validation. We apply the same transformation as previously to pre-process hippocampus volume, ventricles volume and to CSF P-TAU (Box-Cox transform, quantile clipping and min-max rescaling) as well as MMSE and ADAS-Cog13 (only affine-rescaling since scores are already bounded). For GB, we take the Gaussian kernel with σ = 0.24, N<sub>comp</sub> = 4 and g<sub>init</sub> = g<sup>sig</sup> with 2 principals directions for space-shifts, but 3 for DCM. This setting worked best. Nevertheless, we think that it would be interesting to select σ dimension-wise as multidimensional data may not respect the isotropy of a standard Gaussian kernel. In that way, we can choose a smaller σ for features with faster variations and conversely.

**Results** To know whether our results are statistically significant in the absence of other chal-
lengers' errors sample, we derive confidence intervals of DCM and GB's MAE by bootstrapping our
error sample (10000 discounted draws). On figure 1.7, we observe that DCM and GB outperform
the best challenger on ADAS-Cog13 and stay competitive on ventricles volume. This fact shows the
potential of the geometric approach to offer interpretable and flexible models. Regarding the effect
of learning the metric, GB seems to be better than DCM on ADAS-Cog13 and MMSE but not on
ventricles volume. Comparing with the previous experience, it is likely that GB focuses more on
the AD patients profile compared to MCI and controls. Indeed AD patients are representative of
the disease progression whereas MCI and control subjects are seen by the model as slower disease
progressors (delayed in time with a large time-shift τᵢ and slow-paced with a low acceleration factor
αᵢ). It highlights the limits of increasing the flexibility and encourages to combine GB with mixture
models or ensemble methods in the situation of heterogeneity to control the variability.

<footnote n="53">⁵³M. Nguyen, T. He, L. An, et al., "Predicting Alzheimer's disease progression using deep recurrent neural networks," en, *NeuroImage*, vol. 222, p. 117203, Nov. 2020.</footnote>

<footnote n="54">⁵⁴W. Jung, A. W. Mulyadi, and H.-I. Suk, "Unified Modeling of Imputation, Forecasting, and Prediction for AD Progression," en, in *Medical Image Computing and Computer Assisted Intervention – MICCAI 2019*, D. Shen, T. Liu, T. M. Peters, et al., Eds., ser. Lecture Notes in Computer Science, Cham: Springer International Publishing, 2019, pp. 168–176.</footnote>

<footnote n="55">⁵⁵P. J. Moore, T. J. Lyons, J. Gallacher, et al., "Random forest prediction of Alzheimer's disease using pairwise selection from time series data," en, *PLOS ONE*, vol. 14, no. 2, e0211558, Feb. 2019.</footnote>

<footnote n="56">⁵⁶O. Hansson, J. Seibyl, E. Stomrud, et al., "CSF biomarkers of Alzheimer's disease concord with amyloid-PET and predict clinical progression: A study of fully automated immunoassays in BioFINDER and ADNI cohorts," eng,</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 1.7: Boxplots after bootstrapping the prediction errors on the different scores. (For the boxplot: Whiskers=[2.5,97.5] %, Box=[25,50,75]%)</figcaption></figure>

### 1.2.4 Discussion on Geodesic Bending

Riemannian metric learning applied to mixed-effect models enables to reach a sensible trade-off between flexibility and interpretability in disease progression modeling. The proposed approach allows us to disentangle time and space variability while learning the inter-patient variability and the average trajectory from the data. Without data scarcity, it proves to be efficient on homogeneous cohort for improving predictions on imaging and clinical biomarkers and suggests promising perspectives to handle heterogeneous cohorts with ensemble methods and mixture models. One of the main limitations of this work is that we now have some hyper-parameters that require careful tuning. The development of the GB model can be pursued by searching for practical guidelines for selecting hyper-parameters with empirical estimators.

## 1.3 Conclusion

This chapter introduced the Disease Course Mapping model, providing insights into its geometry by drawing a parallel with a simple linear mixed-effect model. We also discussed about the main limitations of the model. One limitation arises from the rigidity imposed by the choice of the metric, which determines the shapes of the possible trajectories. To address this, we introduced Geodesic Bending, an extension of the DCM where the metric is implicitly learned. With this method, we learn a diffeomorphism to compose with the starting geodesic. This results in a deformation of the initial template trajectory in order to match the observed progression patterns. The benefits of this method are sadly overshadowed by a longer optimization process and the cost of hyper-parameter tuning, which was not necessary for the DCM model. Therefore we will stick to the logistic curves DCM in the next chapters of this thesis, which already encompasses our needs for continuous modelling of the markers in neurodegenerative diseases.
# Chapter 2

## Discrete data

Previously we focused on the geometric aspect of the model, whose aim was to better describe the progression of continuous markers on a Riemannian manifold. However features in a medical cohort are varied in nature, and often non-continuous. Symptoms and psychometric scores play a central role in neurodegenerative disease monitoring, yet they are discrete and cannot be modeled as observations on a Riemannian manifold. We will show in this chapter how to extend the Disease Course Mapping model to discrete data. First, we will provide some context on the type of discrete data that we encounter in our longitudinal cohorts. The next section will introduce the binary model and its application to symptom prediction in Parkinson's disease (PD). Then the model will be extended to ordinal data, and the experiments will showcase its potential by analysing a the most widely used clinical score in PD. Finally we will draw a parallel between the DCM extensions to discrete data and the field of Item Response Theory.

Our contributions in this chapter are the two models extending the DCM to binary and ordinal data, and both models' applications led to publications or communications. The application of the binary model was presented in an oral presentation at the Alzheimer's Disease and Parkinson's Disease conference in 2022. The ordinal model was applied to another neurodegenerative disease called ataxia, with the aim to study the Scale for the Assessment and Rating of Ataxia (SARA). This work, performed in collaboration with Paul Moulaire, delivered a clinical paper<fnref n="1" />. Finally the method presented in this chapter has lead to the article<fnref n="2" />. We included most of its content with additional experiments and explanations.

### Contents of the chapter

<table><tr><td>2.1 Context: discrete data, symptoms and disease scales</td><td>53</td></tr><tr><td>2.2 Binary model</td><td>56</td></tr><tr><td>2.2.1 Equation</td><td>56</td></tr></table>

<footnote n="1">¹P. Moulaire, P. E. Poulet, E. Petit, et al., “Temporal Dynamics of the Scale for the Assessment and Rating of Ataxia in Spinocerebellar Ataxias,” en, *Movement Disorders*, vol. 38, no. 1, pp. 35–44, 2023.</footnote>

<footnote n="2">²P.-E. Poulet and S. Durrleman, “Multivariate disease progression modeling with longitudinal ordinal data,” en, *Statistics in Medicine*, vol. n/a, no. n/a,</footnote>
<table><tr><td>2.2.2</td><td>Application to symptom's prediction</td><td>56</td></tr><tr><td>2.3</td><td>Ordinal model</td><td>61</td></tr><tr><td>2.3.1</td><td>Equation</td><td>61</td></tr><tr><td>2.3.2</td><td>Experiments: Simulation study</td><td>62</td></tr><tr><td>2.3.3</td><td>Application to Parkinson's disease data</td><td>66</td></tr><tr><td>2.3.4</td><td>Discussion</td><td>77</td></tr><tr><td>2.4</td><td>Parallel with Item Response Theory</td><td>77</td></tr><tr><td>2.4.1</td><td>Fisher information</td><td>79</td></tr><tr><td>2.5</td><td>Conclusion</td><td>84</td></tr></table>
### 2.1 Context: discrete data, symptoms and disease scales

Today for most neurodegenerative diseases, it is often very expensive to use sensitive early detection methods such as MRI or lumbar puncture, so these techniques are warranted mostly to confirm suspicion of diagnosis. Of course this last statement is likely to change in the next years, as research is currently heavily focused on finding reliable biomarkers for early detection that are also easy to access (mostly blood-based biomarkers). For instance in Alzheimer's disease, phospho-tau blood assays have shown promising results and will soon be implemented in clinics<fnref n="3" />. However for other less prevalent diseases, and even Parkinson's, it is still a long way before automatic early detection based on biomarkers is gold standard. But then, what is the first thing a general practitioner looks at when trying to diagnose a disease? The answer is simple: symptoms. Symptoms are still at the heart of understanding a disease, because they are manifestations of the state of the patient that can be witnessed in everyday life. When talking about Alzheimer's disease, most people immediately think about the loss of memory. For Parkinson's disease, they think about tremor. Symptoms are used to characterize a disease and follow its progression. They are thus often recorded in the longitudinal follow-up of patients, under various formats. Most often they are noted on an intensity scale with a few levels, and sometimes simply written as a binary outcome present/absent.

Neurodegenerative cohorts almost always include clinical assessments of different functions such as cognition, behavior, mood or motor skills. These clinical assessments have been clinically defined specifically for each disease as surrogates for the disease timeline. The cognitive scores tend to be aggregates of subscores corresponding to several aspects of the disease. For instance in Alzheimer's disease we can mention the Alzheimer's Disease Assessment Score-Cognitive subscale<fnref n="4" />, which is a summary of cognitive capabilities regrouped in four categories: language, praxis, memory and orientation. In Parkinson's disease, the MDS-UPDRS cognitive score<fnref n="5" /> is composed of four parts, assessing respectively non-motor symptoms intensity, motor impact in daily life, motor assessment with tests run by the clinician and finally questions related to treatment impact. These scores are built by experts with domain knowledge.

One assessment is composed of several items, e.g. questions or graded tests. However items sometimes only cover one aspect of the disease which may not be shared by the whole population, whereas their aggregation aims at describing a unified and regular pattern of disease progression. For instance in Parkinson's disease one can identify patients with impulse control disorders<fnref n="6" />, but this concerns only 14% of patients according to a study<fnref n="7" />. This information is stored as one of the many items in one of the parts of the MDS-UPDRS scale. As a more general rule, heterogeneity is better assessed through a multivariate analysis of the items before their aggregation.

The global scores aggregate responses to multiple items assessing an array of functional domains.

<footnote n="3">³C. E. Teunissen, I. M. W. Verberk, E. H. Thijssen, et al., "Blood-based biomarkers for Alzheimer's disease: Towards clinical implementation," en, *The Lancet Neurology*, vol. 21, no. 1, pp. 66-77, Jan. 2022.</footnote>

<footnote n="4">⁴J. K. Kueper, M. Speechley, and M. Montero-Odasso, "The Alzheimer's Disease Assessment Scale-Cognitive Subscale (ADAS-Cog): Modifications and Responsiveness in Pre-Dementia Populations. A Narrative Review," *Journal of Alzheimer's Disease*, vol. 63, no. 2, pp. 423-444,</footnote>

<footnote n="5">⁵C. G. Goetz, B. C. Tilley, S. R. Shaftman, et al., "Movement Disorder Society-sponsored revision of the Unified Parkinson's Disease Rating Scale (MDS-UPDRS): Scale presentation and clinimetric testing results: MDS-UPDRS: Clinimetric Assessment," en, *Movement Disorders*, vol. 23, no. 15, pp. 2129-2170, Nov. 2008.</footnote>

<footnote n="6">⁶A. Sharma, V. Goyal, M. Behari, et al., "Impulse control disorders and related behaviours (ICD-RBs) in Parkinson's disease patients: Assessment using "Questionnaire for impulsive-compulsive disorders in Parkinson's disease" (QUIP)," *Annals of Indian Academy of Neurology*, 2015.</footnote>

<footnote n="7">⁷D. Weintraub, J. Koester, M. N. Potenza, et al., "Impulse Control Disorders in Parkinson Disease: A Cross-Sectional Study of 3090 Patients," *Archives of Neurology*, vol. 67, no. 5, pp. 589-595, May 2010.</footnote>
In contrast to it, we aim to enrich disease progression models by breaking down aggregate scores into
subscores or items. In this multivariate analysis, the inter-subject variability is not only described
by changes in the dynamics of progression but also by varying values for each item at a given
disease stage. Disease course mapping is designed to estimate and disentangle these two sources
of variability in a longitudinal data set. However, this approach has been designed for continuous
measurements. The goal of this chapter is to extend this approach to binary or categorical variables.
The study of such data constitutes the field of Item Response Theory (IRT), and we will draw a
parallel with IRT as we proceed to the description of the model.

In IRT, a lot of the concerns of longitudinal models are not taken into account when dealing with non-continuous data. To our knowledge, there is no method in IRT which checks all the following points:

**Disease timeline:** the model should account for individual variability in ages. Indeed patients can have an early or late onset, and we can refer to them as slow or fast progressors. With an adequate time reparametrization, we can map all individuals to a shared disease timeline where it is easier to compare them. As mentioned before, a common disease timeline helps in understanding the standard disease duration and eventually framing disease stages.

**Multivariate:** the method should make use of the multiple observations in the patients' records and jointly model them. Features should not be considered independent of each other, learning the dependencies between them is a key part in understanding disease patterns.

**Inter-patient variability:** disease heterogeneity being a major focus, the model should provide a flexible framework for individual variability. Therefore there should be individual parameters, which are expected to be expressive enough (thus multidimensional) in order to disentangle different disease mechanisms.

In this work we provide a backbone for the longitudinal analysis of binary and ordinal data with the class of mixed-effect models<fnref n="8" />, and more specifically the disease course mapping model<fnref n="9" />. We will first review the literature of disease progression models and the few attempts at modelling non-continuous data. Then we will introduce the DCM model for binary and ordinal data. In the experiments section we showcase the interest of item modelling rather than using aggregated scores in a simple synthetic experiment. The last part of our work is dedicated to applications of the binary and ordinal models to two medical cohorts of Parkinson's disease.

To look at how binary and ordinal data is handled in disease progression models, we need to refer to the main classes of approaches. One such class includes time-to-event modelling, especially with SuStaIn<fnref n="10" /> which is a method mixing time-to-event to define disease stages and subtyping. SuStaIn has recently been extended to deal with ordinal data as well<fnref n="11" />, however this method is taylored for cross-sectional data and does not account for repeated measurements of individuals. The other common framework for longitudinal modelling is the one of mixed-effects models<fnref n="12" />. However such

<footnote n="8">⁸N. M. Laird and J. H. Ware, "Random-effects models for longitudinal data," eng, *Biometrics*, vol. 38, no. 4, pp. 963-974, Dec. 1982.</footnote>

<footnote n="9">⁹J.-B. Schiratti, S. Allassonnière, O. Colliot, et al., "A Bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations," *The Journal of Machine Learning Research*, vol. 18, no. 1, pp. 4840-4872, Jan. 2017.</footnote>

<footnote n="10">¹⁰A. L. Young, R. V. Marinescu, N. P. Oxtoby, et al., "Uncovering the heterogeneity and temporal complexity of neurodegenerative diseases with Subtype and Stage Inference," en, *Nature Communications*, vol. 9, no. 1, pp. 1-16, Oct. 2018.</footnote>

<footnote n="11">¹¹A. L. Young, J. W. Vogel, L. M. Aksman, et al., "Ordinal SuStaIn: Subtype and Stage Inference for Clinical Scores, Visual Ratings, and Other Ordinal Data," *Frontiers in Artificial Intelligence*, vol. 4, 2021.</footnote>

<footnote n="12">¹²N. M. Laird and J. H. Ware, "Random-effects models for longitudinal data," eng, *Biometrics*, vol. 38, no. 4,</footnote>