# Modélisation et analyse stochastique multi-échelle en dynamique des populations : application à la population folliculaire ovarienne

Guillaume Ballif

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

## To cite this version:

Guillaume Ballif. Modélisation et analyse stochastique multi-échelle en dynamique des populations : application à la population folliculaire ovarienne. Analyse numérique [math.NA]. Institut Polytechnique de Paris, 2022. Français. NNT : 2022IPPAX064. tel-04449460

HAL Id: tel-04449460

https://theses.hal.science/tel-04449460

Submitted on 9 Feb 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Modélisation et analyse stochastique multi-échelle en dynamique de population : application à la population folliculaire ovarienne.

Thèse de doctorat de l'Institut Polytechnique de Paris préparée à l'École polytechnique

École doctorale n°574 École doctorale de mathématiques Hadamard (EDMH)
Spécialité de doctorat : Mathématiques appliquées

Thèse présentée et soutenue à Palaiseau, le 27 septembre 2022, par

GUILLAUME BALLIF

## Composition du Jury :

<table><tr><td>Vincent BANSAYE<br>Professeur, École Polytechnique (CMAP)</td><td>Président</td></tr><tr><td>Éva LÖCHERBACH<br>Professeur, Université Paris 1 (SAMM)</td><td>Rapporteur</td></tr><tr><td>Viet-Chi TRAN<br>Professeur, Université Gustave Eiffel (LAMA)</td><td>Rapporteur</td></tr><tr><td>Marc HOFFMANN<br>Professeur, Université Paris-Dauphine (CEREMADE)</td><td>Examineur</td></tr><tr><td>Frédérique CLÉMENT<br>Directeur de recherche, INRIA Saclay - Île de France</td><td>Directeur de thèse</td></tr><tr><td>Romain YVINEC<br>Chargé de recherche, INRAE (PRC)</td><td>Co-directeur de thèse</td></tr></table>
<figure><img src="image_3.png" /></figure>
Remerciements

Je tiens tout d'abord à remercier mes directeurs de thèse Frédérique Clément et Romain Yvinec de m'avoir accompagné durant cette aventure. Merci de m'avoir proposé ce sujet et pour votre aide pour le financement malgré l'océan qui nous séparait. Merci Frédérique pour ta bienveillance, pour ta rigueur scientifique qui m'a fait énormément progresser (surtout pour la rédaction en anglais) et pour ta disponibilité malgré ma tendance à débarquer à l'improviste dans ton bureau. Merci aussi pour ton aide et tes nombreuses explications sur les follicules ovariens et sur la bibliographie disponible. Enfin merci pour cette tradition du tour du lac qui permettait une pause déjeuner bien plus agréable. Merci Romain pour tes connaissance scientifiques, ton recul mathématique qui m'a permis de surmonter les difficultés rencontrées, pour les discussions devant un tableau et pour les nombreuses astuces que tu m'as transmises. Merci aussi pour ta disponibilité et tes nombreuses visites à Saclay. Enfin je voudrais vous remercier de m'avoir laissé partir au Canada pour 4 mois malgré les difficultés sanitaires. Vous m'avez fait confiance pour cette collaboration et je vous en remercie.

I would also like to thank Léa Popovic for her welcome in Montreal. Thank you for all the time you devoted to me during my two visits and for your confidence in our collaboration.

Je souhaite aussi remercier les deux rapporteurs de ce manuscrit, Éva Löcherbach et Viet-Chi Tran, qui m'ont fait l'honneur d'être rapporteurs de mon travail. Je les remercie pour leur lecture attentive et leurs commentaires. Je remercie aussi chaleureusement Vincent Bansaye et Marc Hoffmann pour avoir accepté de faire partie de mon jury.

Durant ces trois années de thèse, j'ai eu la chance d'être entouré par une équipe extraordinaire : l'équipe MDISIM. Merci à l'ensemble de l'équipe pour son accueil, son ambiance et toutes les discussions que j'ai pu avoir. Cela restera dans ma mémoire comme un environnement de travail exceptionnel. Merci aux membres de l'open space Cécile Patte, Jessica Manganotti, Marija Gusseva, Mathieu Barré, Tiphaine Delaunay, Guilia Merlini et Alice Peyraut pour nos nombreuses discussions (sérieuses ou non) et pour votre soutien quotidien. Merci aussi à Jerome Diaz et François Kimmig pour votre aide technique et pratique. Merci à tous les doctorants et stagiaires que j'ai croisés : Anne-Gaëlle Eveno, Louis-Pierre Chaintron, Louis Fostier, Jules Olaye et Mahdi Manoochehrtayebi. Merci à Philippe Moireau, Sébastien Impériale, Martin Genet, Dominique Chapelle et Radek Chabiniok pour nos discussions et votre intérêt pour mon sujet de thèse. Merci à tout ceux qui ont eu le courage de faire le tour du lac avec moi. Nos discussions m'ont permis de faire une vraie pause déjeuner. Merci à Chloé Giraudet, Cecile Patte et Jean-Marc Allain pour nos midi-jeux pendant le confinement !

Je remercie aussi l'ensemble des membres de l'équipe BIOS pour leur accueil à chacune de mes visites à Nouzilly et pour tout le temps qu'ils ont passés à m'écouter parler de mathématiques. Petite mention spéciale pour Frédéric Jean-Alphonse et Thomas Boulo qui sont venus plusieurs fois me chercher à la gare. Merci pour votre aide.
Je souhaite remercier ma famille et mes parents qui m'ont soutenu et accompagné pendant mes études, avec bienveillance et gentillesse. Je vous remercie pour votre soutien et vos encouragements tout au long de mon parcours mathématique.

Enfin, et surtout, je remercie du fond du cœur Amandine. Merci pour ton soutien sans faille, tes encouragements et tes conseils depuis plusieurs années, et plus particulièrement pendant ma thèse. Merci pour ton intérêt pour mon sujet de thèse et pour ton écoute de mes longues explications mathématiques. Tu as été à mes côtés pendant ces trois années, dans les moments difficiles comme dans les moyens plus joyeux, et je n'y serai probablement pas arrivé sans toi.
<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>
# Table des matières

Table des matières i

Table des figures v

Introduction 1

1 Introduction au développement folliculaire ovarien 1

2 Introduction des outils mathématiques 4

2.1 Modèle de populations structurées 4

2.2 Processus stochastique [21] 5

2.3 Processus stochastique multi-échelle [40] 6

2.3.1 Existence de la limite 7

2.3.2 Identification de la limite 7

2.4 Existence et unicité d'une mesure invariante de probabilité 8

2.5 Équation de Poisson 9

3 Plan des chapitres 10

1 Averaging of a stochastic slow-fast model 13

1 Introduction 15

2 Model presentation and main results 16

2.1 Introduction of the model 16

2.2 Hypotheses 16

2.3 Notations and definitions 17

3 Control and ergodicity of the process 19

3.1 Uniform control of the moments of $X^\epsilon$ and $Y^\epsilon$ 19

3.2 Existence and uniqueness of an invariant probability measure 23

4 Proof of the main result 26

4.1 Relative compactness 26

4.2 Identification of the limit 27

4.3 Uniqueness of the limit 29

5 Numerical illustration and biological application 29

5.1 Linear case 30

5.2 Strategy to compute $\pi_x$ and the limit model 31

5.3 Parameter values 31

5.4 Illustration of the convergence 32

5.5 Biological model outputs 32

6 Conclusion 35

7 Appendix 37
<table>
  <tbody>
    <tr>
      <td>7.1</td>
      <td>Generator of the coupling</td>
      <td>37</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>One additional lemma to upper bound B<sup>e</sup>V<sub>p</sub> (lemma 3.2, subsection 3.1)</td>
      <td>37</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Additional Figures</td>
      <td>39</td>
    </tr>
    <tr>
      <td>8</td>
      <td>Simulation numérique de π</td>
      <td>41</td>
    </tr>
    <tr>
      <td>8.1</td>
      <td>Troncature de l'espace d'état</td>
      <td>42</td>
    </tr>
    <tr>
      <td>8.2</td>
      <td>Méthode via la chaîne de Markov</td>
      <td>44</td>
    </tr>
    <tr>
      <td>8.3</td>
      <td>Illustration</td>
      <td>45</td>
    </tr>
    <tr>
      <td>9</td>
      <td>Exploitation biologique du modèle limite</td>
      <td>47</td>
    </tr>
    <tr>
      <td><b>2</b></td>
      <td><b>Modèle à trois échelles de temps</b></td>
      <td><b>51</b></td>
    </tr>
    <tr>
      <td>1</td>
      <td>Motivation</td>
      <td>53</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Modèle à trois échelles de temps</td>
      <td>53</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Changements d'échelle en taille et en temps</td>
      <td>53</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Notations</td>
      <td>54</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Modèle obtenu</td>
      <td>55</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Hypothèses</td>
      <td>56</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Résultats préliminaires</td>
      <td>57</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Contrôle uniforme de Ψ<sup>e</sup></td>
      <td>57</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Comportement en temps long de l'échelle rapide</td>
      <td>58</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Étude de l'échelle intermédiaire</td>
      <td>59</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Modèle limite</td>
      <td>62</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Compacité</td>
      <td>63</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Identification de la limite</td>
      <td>63</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Unicité de la limite</td>
      <td>67</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Discussion du modèle limite</td>
      <td>67</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Étude de l'erreur d'approximation</td>
      <td>68</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Érreur d'approximation pour deux échelles de temps</td>
      <td>68</td>
    </tr>
    <tr>
      <td>6.1.1</td>
      <td>Décomposition de l'erreur</td>
      <td>69</td>
    </tr>
    <tr>
      <td>6.1.2</td>
      <td>Sélection de la fonction H<sup>e</sup></td>
      <td>69</td>
    </tr>
    <tr>
      <td>6.1.3</td>
      <td>Variation quadratique</td>
      <td>71</td>
    </tr>
    <tr>
      <td>6.1.4</td>
      <td>Passage à la limite</td>
      <td>72</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Érreur d'approximation pour trois échelles de temps</td>
      <td>73</td>
    </tr>
    <tr>
      <td>6.2.1</td>
      <td>Décomposition de l'erreur</td>
      <td>74</td>
    </tr>
    <tr>
      <td>6.2.2</td>
      <td>Contrôle de l'erreur issue de l'échelle rapide</td>
      <td>75</td>
    </tr>
    <tr>
      <td>6.2.3</td>
      <td>Contrôle de l'erreur issue de l'échelle intermédiaire</td>
      <td>76</td>
    </tr>
    <tr>
      <td>6.2.4</td>
      <td>Variation quadratique</td>
      <td>78</td>
    </tr>
    <tr>
      <td>6.2.5</td>
      <td>Passage à la limite</td>
      <td>79</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Discussion de la limite obtenue</td>
      <td>80</td>
    </tr>
    <tr>
      <td><b>3</b></td>
      <td><b>Comparison of the model with experimental data on ovarian follicle population</b></td>
      <td><b>83</b></td>
    </tr>
    <tr>
      <td>1</td>
      <td>Introduction</td>
      <td>85</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Experimental datasets on ovarian follicle populations</td>
      <td>86</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Longitudinal data on ovarian follicles numbers as a function of age</td>
      <td>86</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Additional specifications</td>
      <td>87</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Compartmental model of ovarian follicle dynamics with aging</td>
      <td>88</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Model design</td>
      <td>88</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Model analysis : theoretical identifiability</td>
      <td>90</td>
    </tr>
  </tbody>
</table>
<table><tr><td>3.2.1</td><td>Model reparameterization and reduced model</td><td>90</td></tr><tr><td>3.2.2</td><td>Identifiability of the model parameters of the reduced model</td><td>92</td></tr><tr><td>3.2.3</td><td>Identifiability of the initial conditions of the reduced model</td><td>93</td></tr><tr><td>3.2.4</td><td>Identifiability of the initial system</td><td>93</td></tr><tr><td>4</td><td>Parameter estimation strategy</td><td>94</td></tr><tr><td>4.1</td><td>Parameter renormalization and search intervals</td><td>94</td></tr><tr><td>4.2</td><td>Error models for the three datasets</td><td>96</td></tr><tr><td>4.2.1</td><td>Dataset A: whole life data on follicle numbers</td><td>96</td></tr><tr><td>4.2.2</td><td>Error model for point data comparing controlled and uncontrolled follicle activation (Dataset B)</td><td>97</td></tr><tr><td>4.2.3</td><td>Error model for early life data on the contribution of the first-wave follicle pool (Dataset C)</td><td>98</td></tr><tr><td>4.2.4</td><td>Cost function</td><td>98</td></tr><tr><td>4.3</td><td>Parameter estimation results</td><td>99</td></tr><tr><td>4.3.1</td><td>Best fit parameter values</td><td>99</td></tr><tr><td>4.3.2</td><td>Practical identifiability of the parameters</td><td>102</td></tr><tr><td>4.4</td><td>Statistical model selection</td><td>103</td></tr><tr><td>5</td><td>Biological interpretation and model predictions</td><td>104</td></tr><tr><td>5.1</td><td>Insight into the germ cell dynamics</td><td>104</td></tr><tr><td>5.2</td><td>Towards experimental design: an instance of targeted data acquisition counseling</td><td>105</td></tr><tr><td>5.3</td><td>In silico experiments</td><td>106</td></tr><tr><td>6</td><td>Discussion</td><td>108</td></tr><tr><td>7</td><td>Acknowledgements</td><td>111</td></tr><tr><td colspan="2">Discussion et perspectives</td><td>113</td></tr><tr><td colspan="2">Bibliographie</td><td>115</td></tr></table>
<figure><img src="image_6.png" /></figure>
Table des figures

<table>
  <tbody>
    <tr>
      <td>1</td>
      <td>Différents stades de maturation d'un follicule</td>
      <td>2</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Distribution des follicules en fonction de leur diamètre chez trois espèces à différents âges. Figure issue de [63]</td>
      <td>3</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Modèle à d + 1 compartiments avec des taux de migration, de mort et de renouvellement.</td>
      <td>4</td>
    </tr>
    <tr>
      <td>1.1</td>
      <td>Main steps of the limit model simulation</td>
      <td>31</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td>Convergence to the limit model</td>
      <td>32</td>
    </tr>
    <tr>
      <td>1.3</td>
      <td>Comparison between the nonlinear model and three linear formulations fitted on simulated values of x(t) and Y<sub>i</sub>(t)<sub>i=1,2,3</sub>.</td>
      <td>34</td>
    </tr>
    <tr>
      <td>1.4</td>
      <td>Expectation, variance and coefficient of variation of the normalized number of growing follicles as a function of x in the linear and nonlinear cases.</td>
      <td>34</td>
    </tr>
    <tr>
      <td>1.5</td>
      <td>Path from v = (1, 5, 5, 6) to w = (3, 4, 6, 9)</td>
      <td>39</td>
    </tr>
    <tr>
      <td>1.6</td>
      <td>Scheme of the model with three compartments, (X<sup>ε</sup>, Y<sub>1</sub><sup>ε</sup>, Y<sub>2</sub><sup>ε</sup>, Y<sub>3</sub><sup>ε</sup>), used to perform the numerical illustrations in section 5.</td>
      <td>40</td>
    </tr>
    <tr>
      <td>1.7</td>
      <td>Λ<sub>0</sub> (theorem 2.1) as function of x computed with Algorithme 1.</td>
      <td>40</td>
    </tr>
    <tr>
      <td>1.8</td>
      <td>Distribution of the growing follicles at different times</td>
      <td>41</td>
    </tr>
    <tr>
      <td>1.9</td>
      <td>Dynamique du système calculé à l'aide de l'algorithme de troncature</td>
      <td>46</td>
    </tr>
    <tr>
      <td>1.10</td>
      <td>Zoom sur les niveaux N<sub>9</sub>, N<sub>10</sub> et N<sub>11</sub> de la distribution π de la Figure 1.9</td>
      <td>46</td>
    </tr>
    <tr>
      <td>1.11</td>
      <td>Nombres de follicules estimés dans chaque compartiment.</td>
      <td>48</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Follicle numbers as a function of age (Dataset A).</td>
      <td>87</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Additional experimental data.</td>
      <td>88</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Schematic view of the compartmental model (S).</td>
      <td>89</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Logarithmic transformation of follicle numbers in each class as a function of age and spline-based interpolation.</td>
      <td>97</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Fitted follicle numbers in each compartment.</td>
      <td>99</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Fitting results compared with Datasets B and C.</td>
      <td>100</td>
    </tr>
    <tr>
      <td>3.7</td>
      <td>Representative likelihood estimation profiles parameters.</td>
      <td>103</td>
    </tr>
    <tr>
      <td>3.8</td>
      <td>Synthetic data generated with the best fit parameters (Table 3.3).</td>
      <td>105</td>
    </tr>
    <tr>
      <td>3.9</td>
      <td>Estimation performed with the original or augmented datasets.</td>
      <td>107</td>
    </tr>
    <tr>
      <td>3.10</td>
      <td>Likelihood profiles of some parameters for the estimation with synthetic data.</td>
      <td>108</td>
    </tr>
    <tr>
      <td>3.11</td>
      <td>In silico experiment of long-term AMH administration inspired from [38, 37].</td>
      <td>109</td>
    </tr>
  </tbody>
</table>
<figure><img src="image_7.png" /></figure>
# Introduction

Cette thèse porte sur la modélisation mathématique du processus de folliculogenèse tout au long de la vie d'un individu.

## 1 Introduction au développement folliculaire ovarien

Chez les mammifères, l'ovogenèse est le processus aboutissant à la formation et la maturation d'un gamète femelle (appelée ovocyte) permettant la reproduction. Elle est composé de 3 étapes principales : une période de multiplication et différentiation des cellules germinales durant la phase embryonnaire, la formation et la croissance d'un ensemble de cellules (appelé follicule) contenant la cellule germinale et enfin la maturation du follicule jusqu'à l'ovulation (folliculogenèse). La période embryonnaire permet la constitution d'une réserve de follicules quiescents. Chez les grands mammifères (la femme par exemple), cette réserve de follicules est totalement constituée peu de temps après la naissance et va lentement se vider durant l'ensemble de la vie. Chez d'autres mammifères (la souris par exemple), cette réserve se constitue pendant quelques jours juste avant la naissance.

Le développement des follicules ovariens est un long processus de maturation caractérisé par plusieurs stades anatomiques et fonctionnels [56, 49], débutant par l'activation d'un follicule de la réserve et s'achevant soit par la mort au cours de la croissance soit par l'ovulation.

Durant la folliculogenèse, deux phénomènes ont lieu en même temps : l'augmentation de la taille de l'ovocyte (composante germinale) et l'accroissement du nombre de cellules entourant l'ovocyte (composante somatique). Ces cellules s'organisent en couches ce qui permet de définir différents stades de maturation (Figure 1). Peu de follicules atteignent effectivement l'ovulation, la plupart meurt au cours de leur maturation : chez la femme environ 450 follicules participent à l'ovulation sur plus d'un million de follicules primordiaux.

La folliculogénèse peut être séparée en deux phases de croissance du follicule : la croissance basale et la croissance terminale (Figure 1). Comme de nombreux follicules meurent au cours de leur croissance (atrésie) la croissance basale concerne un plus grand nombre de follicules et est plus longue que la croissance terminale. Cette dernière fait de plus intervenir des éléments non ovariens et est fortement liée au cycle ovarien. Dans cette thèse, nous nous concentrerons exclusivement sur l'étude de la croissance basale.

De nombreux signaux et hormones régulent le développement des follicules, certaines étant sécrétées par les follicules eux-mêmes. Les follicules préantraux et les petits follicules antraux (étapes 3 et 4 sur la Figure 1) sécrètent par exemple l'hormone anti-müllerienne (AMH en anglais) qui ralentit l'activation des follicules quiescents [19]. De manière générale, les hormones agissant sur la folliculogenèse ont deux effets complémentaires : le ralentissement de l'activation et de la migration des follicules permettant ainsi d'éviter un épuisement précoce de la réserve,
<figure><img src="image_8.png" /><figcaption>**FIGURE 1** – Différents stades de maturation d'un follicule

Au dessus : Coupes Histologiques de follicules ovariens humains, adapté de [31]
En dessous : Différentes étapes composant la folliculogenèse [49].</figcaption></figure>

et une amplification des taux de mort (principalement dans les stades terminaux) induisant une compétition entre les follicules terminaux et l'ovulation d'un seul follicule (ou d'un petit nombre de follicules chez certains mammifères).

Au cours de la vie, la fonction ovarienne peut être caractérisée à la fois par le nombre de follicules quiescents restants ainsi que par la distribution des follicules selon leur stade de développement [11]. Cette distribution est directement liée aux régulations exercées sur les follicules par les hormones impliquées directement ou indirectement dans la folliculogenèse. La forme de la distribution des follicules est commune à plusieurs espèces malgré les différences d'ordre de grandeur de la réserve de follicules quiescents (Figure 2). Cela motive l'utilisation d'un modèle unique de folliculogenèse pour tous les mammifères permettant à la fois de caractériser le lent épuisement de la réserve et la répartition des follicules entre les différents stades de maturité.

La gestion du nombre et du développement des follicules est crucial pour le potentiel repro-
ductif des femmes. La folliculogenèse étant très sensible à des perturbations d'origine naturelle ou
toxicologique, une modélisation fine de ce processus permettrait de mieux comprendre l'origine et
les conséquences de ces perturbations, voire d'envisager des traitements pour maintenir la fonction
<figure><img src="image_9.png" /></figure>

**FIGURE 2** – Distribution des follicules en fonction de leur diamètre chez trois espèces à différents âges. Figure issue de [63]

ovarienne.

Plusieurs expériences biologiques concernant la population folliculaire ont été effectuées afin d'obtenir des informations qualitatives et quantitatives sur la dynamique folliculaire au cours de la vie. Le dénombrement précis des follicules nécessite d'extraire un ovaire, d'en faire plusieurs coupes puis de compter (souvent à la main) le nombre de follicules présents dans chaque stade de maturation. De telles mesures sont donc très invasives et ne peuvent être faites plusieurs fois sur un même individu. Pour des raisons éthiques, elles ne sont réalisées chez la femme que dans des cas très particuliers. Plusieurs autres espèces (rat, souris, brebis ...) ont été utilisées mais aucunes données de suivi temporel précis de la distribution folliculaire n'existe, excepté chez la souris [34]. Ainsi la plupart des études fournissent des informations partielles sur la folliculogenèse ou bien des informations qualitatives sur la dynamique du développement folliculaire (via des modifications génétiques chez la souris) [19, 66].

De ce fait, la modélisation de la population folliculaire et de son évolution au cours de la vie est cruciale pour mieux comprendre son fonctionnement et pour tenir compte du mieux possible des données disponibles malgré leur hétérogénéité.

Plusieurs modèles ont déjà été développés afin de modéliser une partie ou l'ensemble de la croissance des follicules ovariens chez la souris. Dans [23, 25], les auteurs classifient (à partir des données de [34]) les follicules dans 5 compartiments selon des caractères anatomiques ou physiologiques. Ils considèrent alors un modèle linéaire où chaque follicule est caractérisé par son stade de maturité et peut, soit passer au stade de maturité suivant, soit mourir. Enfin pour tenir compte d'une différence de dynamique entre le début et le reste de la vie, ils introduisent un âge de switch à 20 jours : un jeu de paramètres est estimé sur l'intervalle [0, 20] et un autre sur l'intervalle [20, 600]. Plus récemment, dans [6], les auteurs introduisent deux modèles : le premier, comprenant un nombre fini de compartiments, basé sur des équations différentielles ordinaires ; et le second, considérant la maturation des follicules comme une variable continue, basé sur des équations aux dérivées partielles. Puis tenant compte d'une différence d'échelle de temps au sein de la population
folliculaire, les auteurs conjecturent l'existence d'un modèle limite. Ce modèle limite tient compte
à la fois de la lente exhaustion de la réserve de follicules quiescents ainsi que de la répartition des
follicules au sein des différents stades de maturation.

Dans le chapitre 1, nous prouvons rigoureusement l'existence et l'unicité de ce modèle limite en
utilisant les outils mathématiques introduit dans la section suivante.

2 Introduction des outils mathématiques

2.1 Modèle de populations structurées

Les modèles de populations structurées permettent de suivre l'évolution au cours du temps d'une population. Les individus (cellules, espèces chimiques, follicules...) sont classés selon une variable structurante telle que la taille, le type ou l'âge. Cette variable est a priori continue mais afin de simplifier le modèle, la population peut être structurée en groupes (ou compartiments) correspondant chacun à une valeur différente. Dans toute la suite, nous considérons seulement des modèles de populations structurées avec un nombre fini de compartiments (variable structurante discrète).

<figure><img src="image_10.png" /><figcaption>FIGURE 3 – Modèle à $d + 1$ compartiments avec des taux de migration, de mort et de renouvellement.</figcaption></figure>

La Figure 3 illustre la structure générale d'un modèle à compartiments : la population est divisée en $d+1$ compartiments numérotés de 0 à $d$. Un individu dans le compartiment $k$ peut migrer vers le compartiment suivant (excepté s'il est dans le dernier) ou mourir et ainsi sortir du système. Chacun de ces changements d'état est caractérisé par un taux individuel : $\lambda$ pour la migration et $\mu$ pour la mort. Selon les cas, nous pouvons aussi considérer un taux de renouvellement $r_0$ au sein de la population. Les interactions entre les individus se traduisent par une dépendance de certains de ces taux individuels à l'état de la population.

Dans toute la suite, afin de simplifier les objets mathématiques manipulés, nous faisons deux hypothèses sur la population étudiée :

* les individus sont considérés comme interchangeables au sein d'un compartiment ce qui permet de définir un taux moyen pour chaque compartiment,

* le comportement de la population à l'instant $t$ dépend seulement de l'état de la population à cet instant et est indépendant de l'évolution passée (population sans-mémoire).

Deux points de vue peuvent être adoptés pour interpréter mathématiquement ces modèles : une approche individuelle menant à des processus stochastiques ou une approche moyennée se traduisant par des équations différentielles. Ces modèles ont été énormément appliqués à l'épidémiologie
et à la dynamique cellulaire. Dans le contexte de la folliculogénèse, la principale différence avec ces modèles est l’absence de naissance au sein de la population considérée.

Dans cette partie, nous définissons les différents concepts mathématiques utilisés dans les différents chapitres. La plupart de ces concepts sont reliés aux processus stochastiques et plus particulièrement aux processus stochastiques multi-échelle.

## 2.2 Processus stochastique [21]

Grâce à l'hypothèse de population "sans-mémoire", une première façon d'interpréter mathématiquement un modèle de population structurée fait intervenir les chaînes de Markov à temps continu :

**Definition 2.1**

Un processus stochastique à temps continu $(X(t))_{t \ge 0}$ (ayant pour espace d'état $\mathcal{E}$) est un processus de Markov à temps continu s'il vérifie la propriété de Markov, c'est-à-dire pour tout entier $n \ge 1$,

$$ \mathbb{P}(X(t) = w | X(s) = v, X(t_{n-1}) = v_{n-1}, X(t_{n-2}) = v_{n-2}, \dots, X(t_1) = v_1) = \mathbb{P}(X(t) = w | X(s) = v) \quad (1) $$

avec :

* $0 \le t_1 \le \dots \le t_{n-1} \le s \le t$ est une suite croissante de temps,

* $v_1, v_2, \dots, v_{n-1}, v, w$ sont $n+1$ états parmi l'espace d'état $\mathcal{E}$.

Dans la suite, les processus stochastiques étudiés ont pour espace d'état $(\mathbb{R}^+)^d$ ou $\varepsilon\mathbb{N}^d$ avec $\varepsilon$ un réel strictement positif.

Un processus de Markov à temps continu est totalement défini par son semi-groupe de transition $P = (P_t)_{t \ge 0}$ :

**Definition 2.2**

Soit $f : \mathcal{E} \to \mathbb{R}$ bornée, le semi-groupe $P$ associé au processus de Markov $(X(t))_{t \ge 0}$ est donné par,

$$ P_t f(v) = \mathbb{E}[f(X(t)) | X(0) = v] $$

En pratique, le calcul explicite du semi-groupe $P$ est difficile voire impossible dans certains cas. En effet, le nombre de changement d'états subis par le processus pendant un temps $t$ peut être aussi grand que l'on veut avec une probabilité strictement positive. Afin de pouvoir caractériser le processus indépendamment du temps, l'opérateur suivant est défini :

**Definition 2.3 (Générateur infinitésimal)**

Soit $\mathcal{D}(A) = \{f \text{ bornée, } \lim_{t \to 0} \frac{P_t f - f}{t} \text{ existe pour tout état } x \in \mathcal{E}\}$.

Le générateur $\mathcal{A}$ du processus de Markov à temps continu $(X(t))_{t \ge 0}$ est défini pour tout $f \in \mathcal{A}$.
$$ \mathcal{D}(A), \qquad \mathcal{A}f = \lim_{t \to 0} \frac{P_t f - f}{t} $$

Le générateur infinitésimal ne dépend plus du temps mais seulement des transitions des processus. Ainsi, le processus de Markov à temps continu représentant le modèle à compartiments de la Figure 3 a le générateur suivant, pour toute fonction $g$ bornée à l'infini :

$$ \mathcal{A}g(y) = r_0(y)y_0[g(y + e_0) - g(y)] + \sum_{i=0}^{d} \mu_i(y)y_i[g(y - e_i) - g(y)] \\ + \sum_{i=0}^{d-1} \lambda_i(y)y_i[g(y + e_{i+1} - e_i) - g(y)]. \quad (2) $$

avec $(e_i)_{0 \le i \le d}$ la base canonique de $\mathbb{R}^{d+1}$.

Dans les chapitres 1 et 2, nous introduisons un petit paramètre $\varepsilon$ afin de tenir compte d'une différence d'abondance et de vitesse de transit au sein des compartiments. Le processus ainsi obtenu peut être séparé en plusieurs parties, chacune évoluant à une vitesse temporelle différente. La sous-section suivante détaille la théorie utilisée afin de prouver l'existence et l'unicité d'un modèle limite pour ces processus multi-échelle.

## 2.3 Processus stochastique multi-échelle [40]

Soient $(X^\varepsilon, Y^\varepsilon)_{\varepsilon>0}$ une suite de processus stochastique à valeurs dans $\mathcal{E}_1 \times \mathcal{E}_2$ telle que l'échelle de temps de $Y^\varepsilon$ soit plus rapide que celle de $X^\varepsilon$, c'est-à-dire il existe un opérateur $\mathcal{B}$ tel que pour toute fonction $g$ bornée :

$$ g(Y^\varepsilon(t)) - \int_0^t \beta_\varepsilon \mathcal{B}g(X^\varepsilon(s), Y^\varepsilon(s)) \, ds + \delta_g^\varepsilon(t) \quad (3) $$

est une martingale, avec $\beta_\varepsilon \xrightarrow{\varepsilon \to 0} \infty$ et $\delta_g^\varepsilon/\beta_\varepsilon$ convergant uniformément vers la fonction nulle quand $\varepsilon$ tend vers 0.

De plus, les dynamiques de $X^\varepsilon$ et $Y^\varepsilon$ sont couplées de la façon suivante : il existe un opérateur $\mathcal{C}$ tel que pour toute fonction $f$ bornée :

$$ f(X^\varepsilon(t)) - \int_0^t \mathcal{C}f(X^\varepsilon(s), Y^\varepsilon(s)) \, ds + \alpha_f^\varepsilon(t) \quad (4) $$

est une martingale, avec $\alpha_f^\varepsilon$ convergant uniformément vers la fonction nulle.

L'objectif de la théorie présentée dans cette partie est de prouver (sous des hypothèses raisonnables) que $X^\varepsilon$ converge vers un processus $X$ qui est la solution d'un problème de martingale pour un opérateur obtenu par "moyennisation" de $\mathcal{C}$. Pour cela, comme $Y^\varepsilon$ est plus rapide que la première composante $X^\varepsilon$, la convergence de cette seconde composante n'a pas lieu au sens fonctionnel mais via une mesure. Soit la mesure d'occupation $\Gamma_\varepsilon$ définie pour tout ensemble $B \subset \mathcal{E}_2$ par :
$$ \Gamma_{\varepsilon}([0, t] \times B) := \int_{0}^{t} \mathbb{1}_{B}(Y^{\varepsilon}(s)) \, ds $$

Dans la suite, nous utilisons les notations suivantes :

* $\mathcal{D}_{\varepsilon}[0, \infty[$ : l'ensemble des fonctions à valeurs dans $\mathcal{E}$ qui sont continues à droite et limitées à gauche (càdlàg) sur $[0, \infty[$,

* $\mathcal{L}_m(\mathcal{E})$ : l'ensemble des mesures $\mu$ sur $[0, \infty[ \times \mathcal{E}$ tel que $\mu([0, t[ \times \mathcal{E}) = t$.

### 2.3.1 Existence de la limite

L'existence d'une limite lorsque $\varepsilon$ tend vers 0 découle de la compacité de $(X^{\varepsilon}, \Gamma_{\varepsilon})$ dans $\mathcal{D}_{\varepsilon_1}[0, \infty[ \times \mathcal{L}_m(\mathcal{E}_2)]$. De nombreux critères existent pour la compacité de chacune des composantes dans les espaces appropriés $(\mathcal{D}_{\varepsilon_1}[0, \infty[$ ou $\mathcal{L}_m(\mathcal{E}_2)]$, par exemple un contrôle approprié des moments du processus $(X^{\varepsilon}, Y^{\varepsilon})$ [21].

Dans un second temps, il est nécessaire d'identifier la limite ainsi obtenue.

### 2.3.2 Identification de la limite

Soit $(x, \Gamma)$ un point limite de $(X^{\varepsilon}, \Gamma_{\varepsilon})$. La première étape consiste à identifier la mesure limite $\Gamma$ en utilisant la différence d'échelle de temps au sein du processus $(X^{\varepsilon}, Y^{\varepsilon})$. En divisant par l'équation (3) par $\beta_{\varepsilon}$ et en intégrant le générateur $\mathcal{B}$ contre la mesure $\Gamma_{\varepsilon}$, le processus suivant

$$ \frac{g(Y^{\varepsilon}(t)) + \delta_{g}^{\varepsilon}(t)}{\beta_{\varepsilon}} - \int_{0}^{t} \int_{\mathcal{E}_{2}} \mathcal{B}g(X^{\varepsilon}(s), y) \Gamma_{\varepsilon}(ds \times dy) $$

est une martingale. En faisant tendre $\varepsilon$ vers 0, la mesure limite $\Gamma$ est caractérisée par :

$$ \int_{0}^{t} \int_{\mathcal{E}_{2}} \mathcal{B}g(x(s), y) \Gamma(ds \times dy) = 0 $$

En décomposant la mesure $\Gamma$ afin de séparer la dépendance en temps de la dépendance en $y$ [40], on obtient :

$$ \Gamma(ds \times dy) = \pi_{x(ds)}(dy) \, ds $$

où $\pi_x$ est solution de $\int_{\mathcal{E}_2} \mathcal{B}g(x, y)\pi_x(dy) = 0$.

Cette dernière condition revient à étudier le comportement en temps long de la chaîne de Markov associée au générateur $\mathcal{B}$ pour $x$ fixé, c'est-à-dire l'existence (voire l'unicité) d'une mesure de probabilité invariante. L'ensemble des outils mathématiques intervenant dans cette étude sont présentés dans la Sous-section 2.4.

Enfin, l'identification de la limite $x$ se fait en passant à la limite dans l'équation (4) :

$$ f(x(t)) - \int_{0}^{t} \int_{\mathcal{E}_{2}} \mathcal{C}f(x(s), y) \pi_{x(ds)}(dy) \quad (5) $$

est une martingale.

Le terme intégral apparaissant ci-dessus est la moyennisation du générateur de la première composante $X^{\varepsilon}$ contre la mesure invariante de la seconde composante.

En notant $\bar{C} = \int_{\mathcal{E}_2} \mathcal{C}f(x, y)\pi_x(dy)$ ce nouveau générateur, défini pour une fonction $f$ bornée de $\mathcal{E}_1$ dans $\mathbb{R}$, la limite $x$ est solution de
$$x(t) = x(0) + \int_{0}^{t} \bar{C}f(x(s)) \, ds$$

Pour résumer, un argument de compacité permet d'obtenir l'existence d'une limite pour la séquence de processus $(X^ε, Y^ε)$ lorsque $ε$ tend vers 0. Puis l'identification d'un point limite nécessite l'étude du comportement en temps long d'une chaîne de Markov (pour la partie rapide) et la caractérisation de la moyennisation du générateur de l'échelle lente.

Dans le chapitre 2, l'étude d'un modèle comportant trois échelles de temps utilise les mêmes étapes moyennant quelques adaptations et arguments supplémentaires pour l'identification de la limite.

## 2.4 Existence et unicité d'une mesure invariante de probabilité

Lors de l'étude du comportement en temps long d'un processus de Markov, nous nous intéressons à l'existence puis à l'unicité d'une mesure invariante, c'est-à-dire une mesure $\pi$ telle que

$$\pi P_t = \pi$$

De nombreux critères ont été étudiés afin de prouver l'existence d'une telle mesure $\pi$. Dans la suite, nous nous intéressons à une chaîne de Markov à valeurs dans $\mathcal{E} = \mathbb{N}^d$ (pour d un entier fixé) et nous présenterons le critère de Foster-Lyapunov utilisé dans le chapitre 1.

Pour pouvoir étudier le comportement en temps long d'un processus de Markov, il est nécessaire que celui-ci se comporte bien lorsque le temps devient très grand. Nous supposerons dans la suite que le processus de Markov $(X)$ est non-explosif, c'est-à-dire :

$$\forall t \ge 0, X_t < \infty$$

Définissons tout d'abord le caractère irréductible d'un processus de Markov,

**Definition 2.4**

Le processus $(X_t)$ est irréductible si pour tout état $y, z \in \mathcal{E}$ :

$$\mathbb{E}_y \left( \int_0^\infty \mathbb{1}_{Y_t=z} \, dt \right) > 0.$$

Cette propriété signifie que partant de n'importe quel état, le temps moyen passé par le processus dans chaque état est strictement positif. Un processus non-explosif irréductible est dit Harris recurrent et il admet alors une unique mesure invariante $\pi$ (à une constante près) [2].

Si la mesure est finie, elle peut être renormalisée en une mesure de probabilité. Un processus ayant une unique mesure de probabilité invariante est appelé positif Harris récurrent.

Le critère de Foster-Lyapunov permet de prouver le caractère positif d'un processus Harris récurrent et utilise la notion suivante :

**Definition 2.5**

La fonction $V \ge 0$ est une fonction de Lyapunov pour le générateur $A$ si :
Il existe un ensemble fermé $C \subset \mathcal{E}$ et $a_1, a_2 > 0$ tels que :

$$ \forall y \in \mathcal{E}, \quad AV(y) \leq -a_1[V(y) + 1] + a_2 \mathbb{1}_C(y). $$

Ainsi le théorème suivant permet de conclure sur l'existence et l'unicité d'une mesure de probabilité invariante [48] :

**Théorème 2.1**

Soit $X$ un processus de Markov càdlàg non-explosif. Supposons que $X$ a une fonction de Lyapunov $V$ à valeurs positives avec $C$ un ensemble compact et $V$ bornée sur $C$. Alors le processus $X$ est Harris récurrent positif et

$$ \pi(V) := \sum_{y \in \mathcal{E}} V(y)\pi(y) < \infty $$

Ce théorème ne donne aucune information sur la distribution stationnaire du processus. En particulier, cette distribution dépend des différents paramètres intervenant dans les taux de transition du processus.

## 2.5 Équation de Poisson

Dans la sous-section 2.3, sous les hypothèses appropriées, un processus multi-échelle $(X^\epsilon, Y^\epsilon)$ a pour limite le couple $(x, \pi_x)$ où $x$ est solution de l'équation (5) et $\pi_x$ est la distribution stationnaire d'une chaîne de Markov dépendant de $x$. Afin de prouver l'unicité de ce modèle limite, il est nécessaire de vérifier le caractère lipschitz d'une fonction $S_f(x)$ de la forme suivante, pour $f$ une fonction bornée :

$$ S_f(x) = \langle f, \pi_x \rangle := \sum_{y \in \mathcal{E}} f(y)\pi_x(y) $$

Pour cela, nous avons besoin de caractériser la sensibilité de $\pi_x$ au paramètre $x$ et plus particulièrement d'étudier le **produit scalaire** $\langle f, \pi_x \rangle$. Ainsi il est intéressant d'introduire l'équation de Poisson associée à une fonction $f$ suffisamment régulière par rapport à la fonction de Lyapunov $V$ pour le générateur $A_x$ (dépendant d'un paramètre $x$) :

**Definition 2.6**

Soit une fonction $f$ controlée par la fonction de Lyapunov, c'est-à-dire $|f| \leq c(V + 1)$ pour $c > 0$. L'équation de Poisson associée à $(f, A)$ est l'équation suivante ayant pour inconnue la fonction $h$ :

$$ Ah = \langle f, \pi \rangle - f, \qquad (6) $$

Le théorème suivant traite de l'existence de solution à l'équation de Poisson :

**Théorème 2.2**

Soit $X$ un processus de Markov càdlàg non-explosif. Supposons que $X$ a une fonction de Lyapunov $V$ à valeurs positives avec $C$ un ensemble fermé "petit". Alors le processus $X$ est Harris récurrent positif avec $\pi(V) < \infty$.
De plus, il existe $c < \infty$ tel que pour tout $f \le 1 + V$, l'équation de Poisson associée à $(f, A)$
admet une solution $h$ vérifiant $h \le c(V + 1)$.

Il est à noter que ce théorème est une extension du Théroème 2.1 car il nécessite une hypothèse plus forte pour l'ensemble C :

**Definition 2.7 (Ensemble "petit")**

Soient $a$ une mesure de probabilité sur $\mathbb{R}^+$ et $K_a$ le noyau de transition défini par

$$
K_a = \int_0^t a(s) P_s \, ds
$$

Un ensemble $C$ est dit "petit" s'il existe une mesure non-triviale $\nu_a$ sur $\mathcal{B}(\mathcal{E})$ telle que pour tout $y \in C$ et pour tout $A \in \mathcal{B}(\mathcal{E})$,

$$
K_a(y, A) \geq \nu_a(A)
$$

Le caractère "petit" d'un ensemble C signifie que partant de n'importe quel point y de C, le
processus X a une probabilité strictement positive d'atteindre **n'importe quel état** en un temps fini.

Lors de l'étude de l'erreur entre le modèle initial et son modèle limite, l'équation de Poisson apparaîtra de nouveau dans le calcul de la variance de cette erreur (chapitre 2). Il sera alors nécessaire de caractériser la dépendance de la solution de l'équation au paramètre x ce que nous ferons via l'étude du générateur A_x.

3 Plan des chapitres

Dans cette thèse, nous proposons un modèle général de développement folliculaire ovarien prenant en compte les interactions entre les follicules. Nous introduisons une version stochastique de ce modèle comprenant deux échelles de temps distinctes et prouvons l'existence et l'unicité d'un modèle limite (Chapitre 1). Ce modèle est composé d'une équation différentielle ordinaire couplée avec la distribution stationnaire d'une chaîne de Markov.

D'un point de vue numérique, la suite du chapitre 1 détaille et illustre deux algorithmes spécifiques permettant de simuler efficacement le modèle à deux échelles de temps et plus particulièrement la distribution stationnaire. À l'aide de l'algorithme le plus adapté à notre cas, nous illustrons la convergence vers le modèle limite et certaines caractéristiques du modèle intéressantes du point de vue biologique. Ce chapitre rédigé en anglais est tiré de l'article "Averaging of a stochastic slow-fast model for population dynamics : Application to the development of ovarian follicles." [4].

Dans un second temps, afin de mieux prendre en compte l'accélération progressive de la maturation des follicules, nous ajoutons une troisième échelle de temps au modèle (Chapitre 2). À l'aide des mêmes outils, nous prouvons l'existence et l'unicité d'un modèle limite pour ce nouvel objet mathématique. Il est malgré tout nécessaire de faire des hypothèses plus fortes afin d'obtenir ce résultat. Ce modèle limite est constitué de trois parties couplées entre elles (une pour chaque échelle de temps) : une équation différentielle ordinaire pour l'échelle la plus lente (la réserve de follicules), le point fixe d'une équation différentielle aux dérivées partielles pour l'échelle intermédiaire et la distribution stationnaire d'une chaîne de Markov pour la partie rapide.

Dans la suite du chapitre 2, nous présentons l'étude de l'erreur entre un modèle et sa limite
permettant de quantifier la manière dont chacun de ces modèles convergent vers sa limite. Pour ce faire, il est nécessaire d'étudier l'existence de solutions à différents types d'équations (une équation aux dérivées partielles d'ordre 1 et une équation de Poisson) ainsi que la dépendance de ces solutions à certains paramètres intervenant dans les équations.

Enfin les restrictions du modèle limite ainsi que la difficulté de le simuler dans un contexte biologique réaliste nous ont incités à adapter le modèle initial à deux échelles de temps aux types des données biologiques. Le chapitre 3 rédigé en anglais et soumis au "Journal de Mathematical Biology" présente l'étude de ce modèle. Après s'être assuré de la compatibilité de quatre jeux de données biologiques de nature hétérogène chez la souris, nous étudions l'identifiabilité théorique des paramètres de ce modèle. En nous appuyant sur les résultats de cette identifiabilité, nous avons mis en place une stratégie d'estimation : reparamétrisation du modèle, choix d'un modèle d'errreur et d'une fonction de coût. Les valeurs de paramètres obtenues étant en partie cohérentes avec certaines connaissances biologiques qualitatives, nous avons pu illustrer le potentiel du modèle en reproduisant des expériences in silico.
<figure><img src="image_11.png" /></figure>
# Chapitre 1

## Averaging of a stochastic slow-fast model

This first chapter presents the two-time scale process used to modeling the dynamics of the ovarian follicle population. Taking into account the slow-fast property of the model, we prove the existence and uniqueness of a limit model. Finally, we discuss two methods to simulate it and numerically illustrate the convergence.

Travail réalisé avec Frédérique Clément et Romain Yvinec.

## Sommaire

<table><tr><td>1</td><td>Introduction</td><td>15</td></tr><tr><td>2</td><td>Model presentation and main results</td><td>16</td></tr><tr><td>2.1</td><td>Introduction of the model</td><td>16</td></tr><tr><td>2.2</td><td>Hypotheses</td><td>16</td></tr><tr><td>2.3</td><td>Notations and definitions</td><td>17</td></tr><tr><td>3</td><td>Control and ergodicity of the process</td><td>19</td></tr><tr><td>3.1</td><td>Uniform control of the moments of X<sup>ε</sup> and Y<sup>ε</sup></td><td>19</td></tr><tr><td>3.2</td><td>Existence and uniqueness of an invariant probability measure</td><td>23</td></tr><tr><td>4</td><td>Proof of the main result</td><td>26</td></tr><tr><td>4.1</td><td>Relative compactness</td><td>26</td></tr><tr><td>4.2</td><td>Identification of the limit</td><td>27</td></tr><tr><td>4.3</td><td>Uniqueness of the limit</td><td>29</td></tr><tr><td>5</td><td>Numerical illustration and biological application</td><td>29</td></tr><tr><td>5.1</td><td>Linear case</td><td>30</td></tr><tr><td>5.2</td><td>Strategy to compute π<sub>x</sub> and the limit model</td><td>31</td></tr></table>
<table><tr><td>5.3</td><td>Parameter values</td><td>31</td></tr><tr><td>5.4</td><td>Illustration of the convergence</td><td>32</td></tr><tr><td>5.5</td><td>Biological model outputs</td><td>32</td></tr><tr><td>6</td><td><strong>Conclusion</strong></td><td><strong>35</strong></td></tr><tr><td>7</td><td><strong>Appendix</strong></td><td><strong>37</strong></td></tr><tr><td>7.1</td><td>Generator of the coupling</td><td>37</td></tr><tr><td>7.2</td><td>One additional lemma to upper bound B<sup>e</sup>V<sub>p</sub> (lemma 3.2, subsection 3.1)</td><td>37</td></tr><tr><td>7.3</td><td>Additional Figures</td><td>39</td></tr><tr><td>8</td><td><strong>Simulation numérique de π</strong></td><td><strong>41</strong></td></tr><tr><td>8.1</td><td>Troncature de l'espace d'état</td><td>42</td></tr><tr><td>8.2</td><td>Méthode via la chaîne de Markov</td><td>44</td></tr><tr><td>8.3</td><td>Illustration</td><td>45</td></tr><tr><td>9</td><td><strong>Exploitation biologique du modèle limite</strong></td><td><strong>47</strong></td></tr></table>
# 1 Introduction

We study a birth, migration and death stochastic process modeling the dynamics of a finite population, in which individuals transit unidirectionally across successive compartments. The model is formulated as a continuous-time Markov chain, whose transition matrix involves multiscale effects; the whole (or part of the) population affects the rates of individual birth, migration and death events. The model is characterized by two contrasted time and abundance scales : the first compartment has a high abundance with a slow dynamics, and the remaining compartments have order one abundances with fast dynamics.

The first formulation of the model, as a death and migration process without birth, was introduced in [6] to follow the changes in the population of ovarian follicles throughout lifespan. Ovarian follicles are dynamic endocrine units of the ovaries, enclosing a single germ cell. Each follicle undergoes a sequence of developmental stages, starting with activation from an initially quiescent state, and ending with either ovulation from the final stage or (most of the time) degereneration at any stage. This model is especially suited for mammals, in which a numerous pool of quiescent follicles (the reserve) is constituted once for all in the pre- or peri-natal period. As soon as follicles are formed, and all along life until the pool is exhausted, follicles get slowly activated in an asynchronous way and progress more rapidly through different maturation steps. Hence, at any time, the ovarian function can be characterized by the total number of follicles, and in particular of remaining quiescent follicles, and by the distribution of growing follicles according to their developmental stage. This distribution is further shaped by the control exerted onto the growth (migration) and death rates by reproductive hormones, which emanate (either directly or indirectly) from the follicles themselves (see [11] and references therein). Such controls notably prevent the early exhaustion of the quiescent pool, and are at the source of the cyclic activity of the ovaries. Accordingly, and in contrast to former simpler, linear models [24], the nonlinear expressions for the transition rates include information on follicle numbers in each or all stages. The shape of the follicle distribution is remarkably similar between different mammalian species and with time (except at the very beginning and very end of follicle development, when some compartments are still not filled, or on the contrary, have been exhausted). This distribution both results from and participates in the endocrine homeostasis of the ovaries, and more generally of the reproductive (hypothalamic-pituitary-gonadal) axis. It also ensures that a proper number of mature germ cells is supplied each ovarian cycle for ovulation during all reproductive life.

In [6], taking advantage of the timescale difference between follicle activation and follicle growth, the authors conjectured the existence of a reduced deterministic limit model. In the current study, we derive rigorously the reduction on an extended model including birth events in the reserve. The birth events enable the quiescent pool to be renewed, as encountered in several zoological classes (fishes for instance). We also illustrate numerically the convergence of the full model towards the limit model, as the scaling parameter decreases to zero, and we further explore the behavior of the limit model to highlight the salient features of the nonlinear formulation of the model. The limit model consists of an ordinary differential equation ruling the dynamics of the first compartment, coupled with a quasi-stationary distribution in the remaining compartments, which averages the contribution of the fast component of the Markov chain on the slow one. To derive the limit model, we make use of the generic framework of averaging, suitable to handle singular perturbation problems [40]. Such reduction techniques have found more and more applications in biology, in particular in chemical reaction networks or population dynamics [35]. The main ingredients needed to identify and prove the limit behavior rely on compactness arguments, that we handle here thanks to a Foster-Lyapunov criterion [48], and on coupling arguments [42]. An ad-
ditional crucial step, rarely addressed in similar studies up to our knowledge, is the uniqueness of the solution of the limit model, which we handle by means of the ergodicity theory and Poisson equation, a classical tool in sensitivity analysis [20].

In the linear case (with constant, uncontrolled event rates), direct biologically relevant information can be extracted from the analytical expression of the limit model, as the speed of exhaustion of the first compartment or the shape of the follicle distribution and its sensitivity to parameters. This is not the case in the nonlinear case, for which we settled an appropriate numerical strategy to capture both the dynamics of the slow compartment and the quasi-stationary distribution of the fast compartments. We have combined a simple Euler scheme, to follow the dynamics of the deterministic part, with the nested algorithm [65], to sample efficiently the quasi-stationary distribution of the fast part of the Markov chain. This numerical framework allows us to illustrate the convergence of the initial model towards the limit model, as well as the effects of nonlinear event rates on the model behavior.

The manuscript is organized as follows. The model and the main result are presented in section 2, together with a detailed outline of the main steps for the proof of the main result. Compactness and ergodicity properties through the Foster-Lyapunov criterion are given in section 3. Section 4 is dedicated to the proof of the main result using the averaging tools. In section 5, we illustrate numerically the convergence result, using a nested algorithm to simulate the limit model, and we develop a specific application in reproductive biology. Finally, we discuss our results in section 6.

# 2 Model presentation and main results

## 2.1 Introduction of the model

We consider a population of $X$ quiescent individuals in a slow compartment (the reserve), and a population of individuals distributed into successive fast compartments with increasing maturity $Y = (Y_i)_{i \in [1; d]}$. Individuals may leave the reserve to enter the first fast compartment $Y_1$, at rate $\lambda_0$, or die at rate $\mu_0$. In each compartment $i$, individuals may either progress to the next compartment $i+1$, at rate $\lambda_i$ (except in the last compartment), or die at rate $\mu_i$. In addition, we complete the initial formulation of [6] by adding a birth term, at rate $r_0$ in the slow compartment $X$. After introducing a small parameter $\varepsilon$, rescaling time from $t$ to $t/\varepsilon$ and the reserve size from $X$ to $\varepsilon X$, and multiplying by $\varepsilon$ the rates $\lambda_0$, $\mu_0$ and $r_0$, the model can be fully characterized by the following state transitions :

## 2.2 Hypotheses

To ensure the existence of the process defined in subsection 2.1 at all times, we make the following hypotheses on the rates of events affecting the reserve :

**Hypothesis 2.1 (Birth rate $r_0$ and migration rate $\lambda_0$)**

$$ \exists R_0 > 0, \quad \forall y \in \mathbb{N}^d, r_0(y) \le R_0, $$

$$ \exists B_0 > 0, \quad \forall y \in \mathbb{N}^d, \lambda_0(y) \le B_0. $$

The existence of the stationary probability distribution is strongly linked to the non-accumulation of individuals in a compartment. Accordingly, we assume that
<table><caption>**TABLE 1.1** – Transitions of the process $(X^\epsilon(t), Y^\epsilon(t))_{t \ge 0} \text{in } \mathbb{R}_+ \times \mathbb{N}^d$. We denote by $e = (e_i)_{i \in [1; d]}$ the canonical basis of $\mathbb{N}^d$. To simplify some computations, we set $\lambda_d = 0$.</caption><thead><tr><th></th><th>Transition</th><th>Rate</th></tr></thead><tbody><tr><td>Birth (reserve)</td><td>(X<sup>&epsilon;</sup>, Y<sup>&epsilon;</sup>) &rarr; (X<sup>&epsilon;</sup> + &epsilon;, Y<sup>&epsilon;</sup>)</td><td><sup>r<sub>0</sub>(Y<sup>&epsilon;</sup>)</sup>&frasl;<sub>&epsilon;</sub> X<sup>&epsilon;</sup></td></tr><tr><td>Maturation (reserve)</td><td>(X<sup>&epsilon;</sup>, Y<sup>&epsilon;</sup>) &rarr; (X<sup>&epsilon;</sup> &minus; &epsilon;, Y<sup>&epsilon;</sup> + e<sub>1</sub>)</td><td><sup>&lambda;<sub>0</sub>(Y<sup>&epsilon;</sup>)</sup>&frasl;<sub>&epsilon;</sub> X<sup>&epsilon;</sup></td></tr><tr><td>Death (reserve)</td><td>(X<sup>&epsilon;</sup>, Y<sup>&epsilon;</sup>) &rarr; (X<sup>&epsilon;</sup> &minus; &epsilon;, Y<sup>&epsilon;</sup>)</td><td><sup>&mu;<sub>0</sub>(Y<sup>&epsilon;</sup>)</sup>&frasl;<sub>&epsilon;</sub> X<sup>&epsilon;</sup></td></tr><tr><td>Maturation, i &isin; [1, d &minus; 1]</td><td>(X<sup>&epsilon;</sup>, Y<sup>&epsilon;</sup>) &rarr; (X<sup>&epsilon;</sup>, Y<sup>&epsilon;</sup> &minus; e<sub>i</sub> + e<sub>i+1</sub>)</td><td><sup>&lambda;<sub>i</sub>(Y<sup>&epsilon;</sup>)</sup>&frasl;<sub>&epsilon;</sub> Y<sub>i</sub><sup>&epsilon;</sup></td></tr><tr><td>Death, i &isin; [1, d]</td><td>(X<sup>&epsilon;</sup>, Y<sup>&epsilon;</sup>) &rarr; (X<sup>&epsilon;</sup>, Y<sup>&epsilon;</sup> &minus; e<sub>i</sub>)</td><td><sup>&mu;<sub>i</sub>(Y<sup>&epsilon;</sup>)</sup>&frasl;<sub>&epsilon;</sub> Y<sub>i</sub><sup>&epsilon;</sup></td></tr></tbody></table>

**Hypothesis 2.2 (Migration rates $\lambda_i$)**

$$ \forall i \in [0, d-1], \forall y \in \mathbb{N}^d, \lambda_i(y) > 0, $$

**Hypothesis 2.3 (Non-accumulation)**

$$ \forall i \in [0, d], \alpha_i := \inf_{y \in \mathbb{N}^d} (\lambda_i(y) + \mu_i(y)) > 0. $$

We also make a technical assumption on all event rates :

**Hypothesis 2.4 (Polynomial rates)**

$$ \exists m_0 \in \mathbb{N}, \exists D_0 > 0, \forall y \in \mathbb{N}^d, \left\{ \begin{array}{ll} \forall i \in [1, d-1], & \lambda_i(y) \le D_0 \left[ 1 + \left( \sum_{j=1}^d y_j \right)^{m_0} \right] \\ \\ \forall i \in [0, d], & \mu_i(y) \le D_0 \left[ 1 + \left( \sum_{j=1}^d y_j \right)^{m_0} \right] \end{array} \right. $$

## 2.3 Notations and definitions

We use the following standard notations :

* $C_b(\mathbb{N}^d) = \{f : \mathbb{N}^d \to \mathbb{R} \text{ bounded}\}$,

* $\mathcal{D}_{\mathbb{N}^d}[0, \infty[$: space of $\mathbb{N}^d$-valued functions that are right-continuous and left-limited (càdlàg) on $[0, \infty[$,

* $\mathcal{L}_m(\mathbb{N}^d)$: space of measure $\mu$ on $[0, \infty[ \times \mathbb{N}^d$ such that $\mu([0, t[ \times \mathbb{N}^d) = t$,

* For $x \in \mathbb{R}$, $(x)^+ = \max(0, x)$,
* For $\nu$ a measure on $\mathbb{N}^d$ and $f: \mathbb{N}^d \mapsto \mathbb{R}$ a function integrable against $\nu$, $\langle \nu, f \rangle := \sum_{x \in \mathbb{N}^d} f(x)\nu(x)$.

Also, the following notation will be used repeatedly :

* we define the following sets :

$$ \begin{cases} \mathcal{P} = \{p \in (\mathbb{N}^*)^d, \forall i \in [1, d-1], p_i > p_{i+1}, \} \\ \mathcal{S} = \{w \in \mathbb{N}^d, \forall i \in [1, d-1], w_i \le w_{i+1}\} \end{cases}, $$

* For $Z$ a $d$-dimensional process in $\mathcal{D}_{\mathbb{N}^d}[0, \infty[, $Z_{tot}$ denotes $\sum_{i=1}^d Z_i$.

Let us define rigorously the birth, death and migration process for any $\varepsilon > 0$ :

**Definition 2.1**

Let $(X^\varepsilon(t), Y^\varepsilon(t))_{t \ge 0}$ be the càdlàg $d$-dimensional Markov

process defined in law by the following generator, for any $g: \mathbb{N}^d \times \mathbb{N} \mapsto \mathbb{R}$ bounded,

$$ A^\varepsilon g(x, y) = \frac{r_0(y)}{\varepsilon} x[g(x + \varepsilon, y) - g(x, y)] + \frac{\lambda_0(y)}{\varepsilon} x[g(x - \varepsilon, y + e_1) - g(x, y)] \\ + \frac{\mu_0(y)}{\varepsilon} x[g(x - \varepsilon, y) - g(x, y)] + \sum_{i=1}^d \frac{\mu_i(y)}{\varepsilon} y_i[g(x, y - e_i) - g(x, y)] \\ + \sum_{i=1}^{d-1} \frac{\lambda_i(y)}{\varepsilon} y_i[g(x, y + e_{i+1} - e_i) - g(x, y)]. \quad (1.1) $$

Let $(\mathcal{F}_t^\varepsilon)_{t \ge 0}$ be the canonical filtration associated with $(X^\varepsilon(t), Y^\varepsilon(t))_{t \ge 0}$.

We can now state our main result, the existence and uniqueness of the limit model obtained in the limit $\varepsilon \to 0$ :

**Theorem 2.1**

For any initial condition $(X^\varepsilon(0), Y^\varepsilon(0)) \in \varepsilon\mathbb{N} \times \mathbb{N}^d$ such that, for some deterministic $x^{in} \in \mathbb{R}^+$,

$$ \forall k > 1, \sup_\varepsilon \mathbb{E}(X^\varepsilon(0)^k) < \infty, \sup_\varepsilon \mathbb{E}(Y_{tot}^\varepsilon(0)^k) < \infty, X^\varepsilon(0) \xrightarrow{a.s} x^{in}, $$

the process $(X^\varepsilon, Y^\varepsilon)$ is relatively compact in $\mathcal{D}_{\mathbb{R}^+}[0, \infty[ \times \mathcal{L}_m(\mathbb{N}^d)$, and has a unique limit process $(x, \pi) \in \mathcal{C}^1(\mathbb{R}^+) \times \mathcal{L}_m(\mathbb{N}^d)$ such that :

$$ \left\{ \begin{array}{ll} \frac{dx}{dt}(t) & = \Lambda_0(x(t))x(t), \quad t > 0, \quad \text{with } x(0) = x^{in} \\ \Lambda_0(x(t)) & = \sum_{y \in \mathbb{N}^d} (r_0(y) - \lambda_0(y) - \mu_0(y))\pi_{x(t)}(y) \end{array} \right. $$

where, for $x \in \mathbb{R}^+$, $\pi_x$ verifies

$$ \forall f \in \mathcal{C}_b(\mathbb{N}^d), \langle L_x f, \pi_x \rangle = 0, $$
with $L_x$ the following generator defined for any function $f \in C_b(\mathbb{N}^d)$:

$$ \forall y \in \mathbb{N}^d, L_x f(y) = \lambda_0(y)x[f(y + e_1) - f(y)] + \sum_{i=1}^{d-1} \lambda_i(y)y_i[f(y + e_{i+1} - e_i) - f(y)] \\ + \sum_{i=1}^{d} \mu_i(y)y_i[f(y - e_i) - f(y)]. \quad (1.2) $$

**Remark (Initial condition $x^{in}$):** Theorem 2.1 remains true with $x^{in}$ a random variable with a finite expectation. The proofs remain the same, but are lighter with a deterministic initial condition.

The proof of Theorem 2.1 is divided in two parts, exposed in sections 3 et 4. The main strategy is to use the tools developed in [40], which adapt the standard compactness/identification of the limit scheme to stochastic singular perturbation problems. The compactness of the slow component $X^\epsilon$ is obtained thanks to standard criteria on $D_{\mathbb{R}^+}[0, \infty[$. The fast component $Y^\epsilon$ does not converge in a functional sense, but its occupation measure does, which leads to the study of a (quasi-)stationary problem associated with the generator defined in (1.2). We deploy this machinery in section 4, with compactness criteria in subsection 4.1 and identification of the limit in subsection 4.2. Both suitable moment control on $(X^\epsilon, Y^\epsilon)$ and ergodicity associated with the generator are proved in section 3, using a relevant change of variables. Looking at the cumulative sums associated with $Y^\epsilon$ rather than at the individual numbers in each compartment allows one to exhibit a linear stochastic process that dominates $(X^\epsilon, Y^\epsilon)$ (by coupling methods), and to obtain a uniform moment control on $(X^\epsilon, Y^\epsilon)$ (Proposition 3.1). Similarly, rewriting the generator in terms of cumulative sums allows one to prove the irreducibility, and to find an appropriate Lyapounov function enabling to deduce the existence and uniqueness of $\pi_x$ with minimal assumptions, as well as the local Lipschitz-continuity in $x$ of integrals against $\pi_x$ (Theorem 3.1).

# 3 Control and ergodicity of the process

## 3.1 Uniform control of the moments of $X^\epsilon$ and $Y^\epsilon$

**Proposition 3.1**

Let $(X^\epsilon(0), Y^\epsilon(0)) \in \epsilon\mathbb{N} \times \mathbb{N}^d$ be an initial condition for the process. With the same assumptions as in theorem 2.1, we have:

$$ \forall k > 1, \forall T > 0, \quad \sup_{\epsilon} \mathbb{E} \left( \sup_{t \le T} |X^\epsilon(t)|^k \right) < \infty, $$

$$ \forall k > 1, \quad \sup_{\epsilon} \sup_{t \ge 0} \mathbb{E} \left( |Y_{tot}^{\epsilon}(t)|^k \right) < \infty. $$

To prove Proposition 3.1, we introduce in definition 3.1 a stochastic process $(M^\epsilon, N^\epsilon)$ whose transition rates are linear with respect to the state variables. Writing the fast component in terms of cumulative sums, and coupling $(X^\epsilon, Y^\epsilon)$ with the linear process $(M^\epsilon, N^\epsilon)$, we can show that $(M^\epsilon, N^\epsilon)$ is a majorizing process for $(X^\epsilon, Y^\epsilon)$. The moments of $(M^\epsilon, N^\epsilon)$ can then be controlled thanks to a Lyapounov function (Lemma 3.2).
**Definition 3.1**

Let $(M^\epsilon(t), N^\epsilon(t))_{t \ge 0}$ be the càdlàg $d$-dimensional Markov process defined in law by the following generator $B^\epsilon$, for any function $g \in C_b(\mathbb{R}^+ \times \mathbb{N}^d)$:

$$ B^\epsilon g(m, n) = \frac{R_0}{\epsilon} m[g(m + \epsilon, n) - g(m, n)] + \frac{B_0}{\epsilon} m[g(m, n + e_1) - g(m, n)] \\ + \frac{\alpha_0}{\epsilon} m[g(m - \epsilon, n + e_1) - g(m, n)] + \sum_{i=1}^{d} \frac{\alpha_i}{\epsilon} n_i [g(m, n + e_{i+1} - e_i) - g(m, n)]. \quad (1.3) $$

We introduce the following change of variables for $1 \le i \le d$,

$$ V_i^\epsilon(t) = \sum_{j=1}^{i} Y_j^\epsilon(t), \quad W_i^\epsilon(t) = \sum_{j=1}^{i} N_j^\epsilon(t). \quad (1.4) $$

It is clear that $V^\epsilon, W^\epsilon \in S$. With a slight abuse of notation, the generators $A^\epsilon$ (1.1) and $B^\epsilon$ (1.3) become for any $v, w \in S$ (with $v_0 = w_0 = 0$):

$$ \begin{aligned} A^\epsilon g(x, v) = & \frac{r_0(v)}{\epsilon} x[g(x + \epsilon, v) - g(x, v)] + \frac{\lambda_0(v)}{\epsilon} x[g(x - \epsilon, v + \sum_{j=1}^d e_j) - g(x, v)] \\ & + \frac{\mu_0(v)}{\epsilon} x[g(x - \epsilon, v) - g(x, v)] \\ & + \sum_{i=1}^d \frac{\mu_i(v)}{\epsilon} (v_i - v_{i-1}) [g(x, v - \sum_{j=i}^d e_j) - g(x, v)] \\ & + \sum_{i=1}^{d-1} \frac{\lambda_i(v)}{\epsilon} (v_i - v_{i-1}) [g(x, v - e_i) - g(x, v)], \end{aligned} \quad (1.5) $$

$$ \begin{aligned} B^\epsilon g(m, w) = & \frac{R_0}{\epsilon} m[g(m + \epsilon, w) - g(m, w)] + \frac{B_0}{\epsilon} m[g(m, w + \sum_{i=1}^d e_i) - g(m, w)] \\ & + \frac{\alpha_0}{\epsilon} m[g(m - \epsilon, w + \sum_{i=1}^d e_i) - g(m, w)] \\ & + \sum_{j=1}^d \frac{\alpha_i}{\epsilon} (w_i - w_{i-1}) [g(m, w - e_i) - g(m, w)]. \end{aligned} \quad (1.6) $$

We prove in the next lemma that $(X^\epsilon, V^\epsilon)$ is dominated by $(M^\epsilon, W^\epsilon)$.

**Lemma 3.1**

For any $(M^\epsilon(0), N^\epsilon(0))$ such that $X^\epsilon(0) \le M^\epsilon(0)$ and $Y^\epsilon(0) \le N^\epsilon(0)$, we have, for any $t > 0$,

$$ \left\{ \begin{array}{l} X^\epsilon(t) \le M^\epsilon(t) \text{ stochastically,} \\ V^\epsilon(t) \le W^\epsilon(t) \text{ stochastically componentwise.} \end{array} \right. $$
**Proof:** We introduce a coupling between processes $(X^\epsilon, V^\epsilon)$ and $(M^\epsilon, W^\epsilon)$ with the following events driven by the slow components $(X^\epsilon, M^\epsilon)$:

$$ \left\{
\begin{array}{ll}
(x, v, m, w) \rightarrow (x + \epsilon, v, m + \epsilon, w) & \text{at rate } \frac{r_0(v)}{\epsilon}x, \\
(x, v, m, w) \rightarrow (x, v, m + \epsilon, w) & \text{at rate } \left(\frac{R_0}{\epsilon}m - \frac{r_0(v)}{\epsilon}x\right)^+, \\
(x, v, m, w) \rightarrow (x - \epsilon, v + \sum_{j=1}^{d} e_j, m, w + \sum_{j=1}^{d} e_j) & \text{at rate } \frac{\lambda_0(v)}{\epsilon}x, \\
(x, v, m, w) \rightarrow (x, v, m, w + \sum_{j=1}^{d} e_j) & \text{at rate } \left(\frac{B_0}{\epsilon}m - \frac{\lambda_0(v)}{\epsilon}x\right)^+, \\
(x, v, m, w) \rightarrow (x - \epsilon, v, m, w) & \text{at rate } \frac{1}{\epsilon}(\mu_0(v)x - \alpha_0 \min(x, m)), \\
(x, v, m, w) \rightarrow (x - \epsilon, v, m - \epsilon, w + \sum_{j=1}^{d} e_j) & \text{at rate } \frac{\alpha_0}{\epsilon}\min(x, m), \\
(x, v, m, w) \xrightarrow{7^{\text{th}}} (x, v, m - \epsilon, w + \sum_{j=1}^{d} e_j) & \text{at rate } \frac{\alpha_0}{\epsilon}(m - \min(x, m)),
\end{array}
\right. $$

and the following events driven by the fast components $(V^\epsilon, W^\epsilon)$, for any $i \in [1, d]$:

$$ \left\{
\begin{array}{ll}
(x, v, m, w) \rightarrow (x, v - e_i, m, w - e_i) & \text{at rate } \frac{1}{\epsilon} \frac{\alpha_i u_i(v, w) \lambda_i(v)}{\lambda_i(v) + \mu_i(v)}, \\
(x, v, m, w) \rightarrow (x, v - \sum_{j=i}^{d} e_j, m, w - e_i) & \text{at rate } \frac{1}{\epsilon} \frac{\alpha_i u_i(v, w) \mu_i(v)}{\lambda_i(v) + \mu_i(v)}, \\
(x, v, m, w) \xrightarrow{10^{\text{th}}} (x, v, m, w - e_i) & \text{at rate } \frac{\alpha_i}{\epsilon}(w_i - w_{i-1} - u_i(v, w)), \\
(x, v, m, w) \rightarrow (x, v - \sum_{j=i}^{d} e_j, m, w) & \text{at rate } \frac{1}{\epsilon} \left[ \mu_i(v)(v_i - v_{i-1}) - \frac{\alpha_i u_i(v, w) \mu_i(v)}{\lambda_i(v) + \mu_i(v)} \right], \\
(x, v, m, w) \rightarrow (x, v - e_i, m, w) & \text{at rate } \frac{1}{\epsilon} \left[ \lambda_i(v)(v_i - v_{i-1}) - \frac{\alpha_i u_i(v, w) \lambda_i(v)}{\lambda_i(v) + \mu_i(v)} \right],
\end{array}
\right. $$

with $u_i(v, w) = \min(v_i - v_{i-1}, w_i - w_{i-1})$.

The marginal $(X^\epsilon, V^\epsilon)$ has for generator $A^\epsilon$ (1.5) and the marginal $(M^\epsilon, W^\epsilon)$ has for generator $B^\epsilon$ (1.6), as long as $X^\epsilon \le M^\epsilon$ (see subsection 7.1). Furthermore, each event preserves the order $X^\epsilon \le M^\epsilon$ except the seventh event. Since $X^\epsilon, M^\epsilon \in \epsilon\mathbb{N}$ and the jump increment of the seventh event is $-\epsilon$, the order between $X^\epsilon$ and $M^\epsilon$ can be reversed only when $X^\epsilon = M^\epsilon$, in which case the rate is zero. Thus, the order $X^\epsilon \le M^\epsilon$ is preserved.

In the same way, each event preserves the (componentwise) order $V^\epsilon \le W^\epsilon$ except the tenth event. For any $i \in [1, d]$, since $V_i^\epsilon, W_i^\epsilon \in \mathbb{N}^d$ and the jump increment of the tenth event is $-1$, the order between $V_i^\epsilon$ and $W_i^\epsilon$ can be reversed only when $V_i^\epsilon = W_i^\epsilon$. However, in any state such that

$$ V_i^\epsilon = W_i^\epsilon, \text{ and } \forall j \neq i, V_j^\epsilon \le W_j^\epsilon, $$

the rate of the tenth event is equal to 0, so that the order $V^\epsilon \le W^\epsilon$ is preserved, which achieves the proof.

We now prove that the moments of $(M^\epsilon, W^\epsilon)$ are uniformly bounded.
**Lemma 3.2**

For any sequence of initial condition $(M^\epsilon(0), W^\epsilon(0)) \in \epsilon\mathbb{N} \times \mathbb{N}^d$ such that

$$ \forall k > 1, \sup_{\epsilon} \mathbb{E}(M^{\epsilon}(0)^k) < \infty, \sup_{\epsilon} \mathbb{E}(W_{tot}^{\epsilon}(0)^k) < \infty, $$

we have, for the slow component,

$$ \forall k > 1, \forall T > 0, \sup_{\epsilon} \mathbb{E} \left( \sup_{t \le T} |M^{\epsilon}(t)|^k \right) < \infty, \quad (1.7) $$

and, for the fast component,

$$ \forall k > 1, \sup_{\epsilon} \sup_{t \ge 0} \mathbb{E} \left( |W_{tot}^{\epsilon}(t)|^k \right) < \infty. \quad (1.8) $$

**Proof:** $M^\epsilon$ is a linear birth and death process with birth (resp. death) rate $\frac{R_0}{\epsilon}m$ (resp. $\frac{\alpha_0}{\epsilon}m$). Since the jumps increment $M^\epsilon$ by $\pm\epsilon$, we obtain (1.7) by a straightforward adaptation of [46, Proposition 2.7]. Let $p \in \mathcal{P}$. We will prove that function $V_p$, defined for any $w \in \mathcal{S}$ by $V_p(w) = \sum_{i=1}^d w_i^{p_i}$, is a Lyapunov function for $W^\epsilon$, that is, there exists a pair $a_1, a_2 > 0$ and a compact set $C \subset \mathbb{N}^d$ such that:

$$ B^{\epsilon} V_{p}(m, w) \leq \frac{a_2}{\epsilon} \mathbb{1}_{C}(w) - \frac{a_1}{\epsilon} V_{p}(w). \quad (1.9) $$

From (1.6) and using lemma 7.1, we obtain, for any $w \in \mathcal{S}$:

$$ \begin{aligned} B^{\epsilon} V_{p}(m, w) &= \frac{B_0 + \alpha_0}{\epsilon} m \sum_{i=1}^{d} \left[ (w_i + 1)^{p_i} - w_i^{p_i} \right] \\ &\quad + \sum_{i=1}^{d} \frac{\alpha_i}{\epsilon} (w_i - w_{i-1}) \left[ (w_i - 1)^{p_i} - w_i^{p_i} \right] \\ &\leq \frac{B_0 + \alpha_0}{\epsilon} m \sum_{i=1}^{d} 2^{p_i} (w_i^{p_i-1} + 1) - \sum_{i=1}^{d} \frac{\alpha_i}{\epsilon} (w_i - w_{i-1}) w_i^{p_i-1} \\ &\leq \frac{F(w)}{\epsilon} - \frac{\alpha}{2\epsilon} V_p(w), \end{aligned} \quad (1.10) $$

where $\alpha = \min_{1 \le i \le d-1} \alpha_i > 0$ from hypothesis 2.3, and $F$ is given by, for $w \in \mathcal{S}$,

$$ F(w) = (B_0 + \alpha_0)m \sum_{i=1}^{d} 2^{p_i} (w_i^{p_i-1} + 1) + \alpha \sum_{i=1}^{d} w_{i-1} w_i^{p_i-1} - \frac{\alpha}{2} V_p(w). $$

To prove (1.9), we observe that, for any $i$, $w_i \le V_p(w)^{1/p_i}$, so that, for any $p \in \mathcal{P}$,

$$ \lim_{\|w\| \to \infty} \frac{F(w)}{V_p(w)} = -\frac{\alpha}{2}. $$

Hence, there is a compact set $C$ such that $F$ is negative on $\mathbb{N}^d \setminus C$. In consequence, $F$ reaches a finite maximal value on $\mathbb{N}^d$, leading to (1.9) with $a_1 = \alpha/2$ and $a_2 = \max_{w \in \mathbb{N}^d} F(w) < \infty$. We now deduce
from (1.9) that

$$ \begin{aligned} \mathbb{E}[V_p(W^\epsilon(t))] &\le \exp\left(-\frac{a_1}{\epsilon}t\right) \mathbb{E}[V_p(W^\epsilon(0))] + \frac{a_2}{a_1} \left(1 - \exp\left(-\frac{a_1}{\epsilon}t\right)\right) \\ &\le \mathbb{E}[V_p(W^\epsilon(0))] + \frac{a_2}{a_1}. \end{aligned} \quad (1.11) $$

Since (1.11) is verified for all $p \in \mathcal{P}$, all $\epsilon > 0$ and all $t \ge 0$, we finally obtain (1.8). ■

**Proof of Proposition 3.1:** We apply lemma 3.1 with lemma 3.2, and use the inequality $Y_{tot}^\epsilon \le V_{tot}^\epsilon$, to obtain Proposition 3.1. ■

## 3.2 Existence and uniqueness of an invariant probability measure

We introduce the process $Y$ which represents the fast part of $(X^\epsilon, Y^\epsilon)$, for a frozen value of $X^\epsilon$. In all this subsection, let $x > 0$ be a fixed positive number.

**Definition 3.2**

Let $Y$ be the càdlàg $d$-dimensional Markov process defined in law by the following generator, for any $f: \mathbb{N}^d \times \mathbb{N} \to \mathbb{R}$ bounded,

$$ \forall y \in \mathbb{N}^d, L_x f(y) = \lambda_0(y)x[f(y + e_1) - f(y)] + \sum_{i=1}^{d-1} \lambda_i(y)y_i[f(y + e_{i+1} - e_i) - f(y)] \\ + \sum_{i=1}^{d} \mu_i(y)y_i[f(y - e_i) - f(y)]. $$

The next theorem states the main result on the ergodicity property of $Y$.

**Theorem 3.1**

$Y$ is exponentially-ergodic and has a unique invariant probability measure $\pi_x$. Furthermore, function $S_f$ defined by

$$ \forall x \in \mathbb{R}^+, S_f(x) = \langle f, \pi_x \rangle, $$

is locally Lipschitz continuous on $\mathbb{R}^+$, if, for some $p \in \mathcal{P}$ and $c > 0$, $|f| \le c(V_p + 1)$. ■

We recall that

$$ \forall p \in \mathcal{P}, V_p(y) = \sum_{i=1}^{d} \left( \sum_{k=1}^{i} y_k \right)^{p_i}. $$

To prove theorem 3.1, we apply the criteria of [48]. Namely, we show below that $Y$ is irreducible (lemma 3.3) by finding a path joining any pair of points in the state space, and we next show that function $V_p$ is a Lyapounov function for $Y$ (lemma 3.4). Similarly to subsection 3.1, the main idea is to look at cumulative sums instead of individual numbers in each compartment.

**Lemma 3.3**

For any $y, z \in \mathbb{N}^d$:

$$ \mathbb{E}_y \left( \int_0^\infty \mathbb{1}_{Y_t=z} dt \right) > 0. $$
**Proof:** Let $y, z \in \mathbb{N}^d$. We now construct a path $C_y^z$ from $y$ to $z$ that involves $N$ events, such that the set

$$A = \{ \text{choose } C_x^y \text{ among paths of N events} \mid \text{at least N events occur} \},$$

has a positive probability. We further show that the time needed to perform at least $N$ events is finite almost surely, leading to the positive lower bound

$$\mathbb{E}_y \left[ \int_0^\infty 1_{Y_t=z} dt \right] \ge \left( \sum_{i=1}^d (\lambda_i(z) + \mu_i(z)) \right)^{-1} \mathbb{P}(A) > 0,$$

since, when $Y_t$ reaches state $z$, it stays in $z$ for a mean time $(\sum_{i=1}^d \lambda_i(z) + \mu_i(z))^{-1} > 0$. To simplify the construction of $C_y^z$, we use the same change of variables as in (1.4), with $v_i = \sum_{j=1}^i y_j$ and $w_i = \sum_{j=1}^i z_j$ (see Figure 1.5):

(1) **Total number compensation step.** Starting from $v$, perform

$$M = \max (0, w_1 - v_1, w_2 - v_2, \dots, w_d - v_d) \text{ times the immigration event of an individual into the first compartment, leading to state}$$

$$v^{(1)} = (v_1 + M, v_2 + M, \dots, v_d + M).$$

These immigration events occur at rate $\lambda_0(\cdot)x > 0$ (hypothesis 2.2). Note that, by definition of $M$, $M \ge w_k - v_k$ for any $k$, and $v^{(1)} \ge w$ componentwise.

(2) **Compartment-balancing step.** For $k$ from 1 to $d-1$, starting from $v^{(k)} = (w_1, \dots, w_{k-1}, v_k + M, \dots, v_d + M)$,

* If $M = w_k - v_k$, then $v^{(k)} = v^{(k+1)}$ and the state is left unchanged.

* If $M > w_k - v_k$, perform $M - (w_k - v_k)$ times the migration event from the $k^{th}$ compartment into the next one to reach state $v^{(k+1)}$. For $i = M + v_k, \dots, w_k + 1$, each migration event occurs at rate $\lambda_k(\cdot)(i - w_{k-1})$. Since $w_k \ge w_{k-1}$ and given hypothesis 2.2, all above rates are strictly positive.

(3) **Surplus-removing step.** From state $v^{(d)} = (w_1, \dots, w_{d-1}, v_d + M)$,

* If $M = w_d - v_d$, then $v^{(d)} = w$ and the path has reached $z$.

* if $M > w_d - v_d$, perform $M - (w_d - v_d)$ times the death event in the last compartment to reach $w$. For $i = M + v_d, \dots, w_d + 1$, each event occurs at rate $\mu_d(\cdot)(i - w_{d-1})$ and are strictly positive (hypothesis 2.3).

The number $N$ of events in $C_y^z$ is clearly finite, and can be bounded by $(d+1)M + \sum v_k$. Since we create no more than $M$ individuals along $C_y^z$, the states stay in a compact of $\mathbb{N}^d$, and the total event rate is upper bounded. As the rate of each event occurring along $C_x^y$ is strictly positive, we have $\mathbb{P}(A) > 0$. Finally, as $Y$ has no absorption point (see hypothesis 2.3), the time needed to perform at least $N$ events is clearly finite almost surely.

We now turn to the Lyapounov property.

**Lemma 3.4**

For $\forall p \in \mathcal{P}$, $V_p$ is a Lyapunov function for the generator $L_x$, that is : There exists a closed set
$C \subset \mathbb{N}^d$ and $a_1, a_2 > 0$ such that :

$$ \forall y \in \mathbb{N}^d, L_x V_p(y) \leq -a_1 [V_p(y) + 1] + a_2 1_C(y). $$

**Proof:** Once again, we apply the change of variables (1.4) on the generator $L_x$ (with a slight abuse of notation) :

$$ \begin{aligned} L_x V_p(w) ={}& \lambda_0(w)x \left[V_p(w + \sum_{i=1}^d e_i) - V_p(w)\right] \\ & + \sum_{i=1}^{d-1} \lambda_i(w)(w_i - w_{i-1}) \left[V_p(w - e_i) - V_p(w)\right] \\ & + \sum_{i=1}^d \mu_i(y)(w_i - w_{i-1}) \left[V_p(w - \sum_{j=i}^d e_j) - V_p(w)\right]. \end{aligned} $$

For any $i \in [1, d]$, $V_p(w - \sum_{j=i}^d e_j) \leq V_p(w - e_i) \leq V_p(w)$, and, given hypotheses 2.1 et 2.3, we have, for any $w \in S$ :

$$ L_x V_p(w) \leq B_0 x \sum_{i=0}^{d} \left[ (w_i + 1)^{p_i} - w_i^{p_i} \right] + \sum_{i=1}^{d} \alpha_i (w_i - w_{i-1}) \left[ (w_i - 1)^{p_i} - w_i^{p_i} \right]. $$

The same algebraic manipulations as in (1.10) lead to lemma 3.4.

**Proof of Theorem 3.1:** Combining lemmas 3.3 et 3.4, with [48, Theorem 4.4], we obtain that $Y$ has a unique invariant probability measure $\pi_x$, and, moreover, for any $p \in \mathcal{P}$,

$$ \langle V_p, \pi_x \rangle < \infty. \tag{1.12} $$

To show that $S_f$ is locally Lipschitz continuous, we first prove that function $S_{V_p}$ is locally bounded, for any $p \in \mathcal{P}$, by a coupling approach. Let $V$ be the cumulative sums associated with $Y$ (see (1.4)), and let $x \in (0, K)$ with $K < \infty$. By analogy with (1.6), let $W$ be the càdlàg d-dimensional Markov process defined in law by the following generator :

$$ Bg(w) = B_0 K \left[ g(w + \sum_{i=1}^{d} e_i) - g(w) \right] + \sum_{j=1}^{d} \alpha_i (w_i - w_{i-1}) \left[ g(w - e_i) - g(w) \right]. $$

We take $W(0) \geq V(0)$. Since $\lambda_0(\cdot)x < B_0 K$, adapting the coupling used to prove lemma 3.1, we have that, for any $t > 0$, for any $i \in [1, d]$, $V_i(t) \leq W_i(t)$ stochastically. Moreover $W$ is also clearly exponentially ergodic, and thanks to Lyapounov function $V_p$, its stationary measure $\nu_K$ satisfies, for any $x < K$, any $p \in \mathcal{P}$,

$$ \langle V_p, \pi_x \rangle \leq \langle V_p, \nu_K \rangle < \infty. $$

Let $f$ be a function such that $|f| \leq c(V_p + 1)$ for some $p \in \mathcal{P}$ and $c > 0$. It remains to show that function $S_f$ is locally Lipschitz continuous on $\mathbb{R}^+$. We introduce the Poisson equation associated with $(f, L_x)$ :

$$ L_x h = \langle f, \pi_x \rangle - f, \tag{1.13} $$
and, for any function $h: \mathbb{N}^d \to \mathbb{R}$, the $(V_p + 1)$-norm

$$ \|h\|_{V+1} := \sup_{x \in \mathbb{N}^d} \frac{|h(x)|}{V_p(x) + 1}. $$

Thanks to the Lyapounov property lemma 3.4, we can apply [29, Theorem 3.2] to obtain that, for any $x > 0$, there exists $\alpha_x < \infty$ such that (1.13) admits a solution, $g_x$, with $|g_x| \le \alpha_x(V_p + 1)$. From (1.12), for any $x_1, x_2 \in \mathbb{R}^+$, $\pi_{x_2}(g_{x_1}) < \infty$, and we have :

$$ \begin{aligned} \langle f, \pi_{x_2} - \pi_{x_1} \rangle &= \langle f, \pi_{x_2} \rangle - \langle f, \pi_{x_1} \rangle \langle 1, \pi_{x_2} \rangle = \langle f - \langle f, \pi_{x_1} \rangle 1, \pi_{x_2} \rangle \\ &= \langle -L_{x_1} g_{x_1}, \pi_{x_2} \rangle = \langle -L_{x_1} g_{x_1} + L_{x_2} g_{x_1}, \pi_{x_2} \rangle. \end{aligned} $$

From (1.2), for any $y \in \mathbb{N}^d$,

$$ -L_{x_1}g_{x_1}(y) + L_{x_2}g_{x_1}(y) = \lambda_0(y)(x_2 - x_1)[g_{x_1}(y + e_1) - g_{x_1}(y)]. $$

Since $g_{x_1}$ has a finite $(V+1)$-norm, it is also true for $\tilde{g}_{x_1} : y \mapsto g_{x_1}(y + e_1)$, hence $\pi_{x_2}(\tilde{g}_{x_1}) < \infty$, and we have

$$ \langle f, \pi_{x_2} - \pi_{x_1} \rangle = (x_2 - x_1)\langle (\tilde{g}_{x_1} - g_{x_1})\lambda_0, \pi_{x_2} \rangle. $$

Let K be a neighborhood of $x_1$, we obtain that, for any $x_2 \in K$,

$$ \left| \langle f, \pi_{x_2} - \pi_{x_1} \rangle \right| \le |x_2 - x_1| \max (\|\tilde{g}_{x_1}\|_{V_p+1}, \|g_{x_1}\|_{V_p+1}) \times 2B_0 \sup_{y \in K} \langle V_p, \pi_y \rangle, $$

which concludes the proof.

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

# 4 Proof of the main result

In this section, we will prove theorem 2.1 in three steps : relative compactness of $(X^\epsilon, Y^\epsilon)$ in an appropriate topology, identification of a limiting equation for any adherence value, and uniqueness of the limit. The first two steps rely on the tools developed in [40] with the estimates provided by Proposition 3.1. The last step relies on theorem 3.1.

## 4.1 Relative compactness

We start by the relative compactness of $X^\epsilon$ in $\mathcal{D}_{\mathbb{R}^+}[0, \infty[$. Let $A : C_b^2(\mathbb{R}^+) \to C_b^1(\mathbb{R}^+ \times \mathbb{N}^d)$ be the infinitesimal generator defined for any $x \in \mathbb{R}^+$ and $y \in \mathbb{N}^d$ by :

$$ Af(x, y) = (r_0(y) - \lambda_0(y) - \mu_0(y))xf'(x). $$

Then, using definition 2.1 with $g$ depending only on the first variable, i.e. $\forall y \in \mathbb{N}^d$, $g(\cdot, y) = f(\cdot)$, we obtain that

$$ M_f^\epsilon(t) = f(X^\epsilon(t)) - f(X^\epsilon(0)) - \int_0^t Af(X^\epsilon(s), Y^\epsilon(s)) ds + R_f^\epsilon(t), \quad (1.14) $$

is a $\mathcal{F}_t^\epsilon$-martingale, where $R_f^\epsilon(t) = \int_0^t [Af - A^\epsilon f](X^\epsilon(s), Y^\epsilon(s)) ds$. Applying [21, Theorems 3.9.1 and 3.9.4], $X^\epsilon$ is relatively compact in $\mathcal{D}_{\mathbb{R}^+}[0, \infty[$ if, for any $f \in C_b^2(\mathbb{R}^+)$:
• $\mathbb{E}(\sup_{t \le T} |R_f^\varepsilon(t)|) \xrightarrow{\varepsilon \to 0} 0,$

• $\exists q > 1$ such that $\sup_{\varepsilon} \mathbb{E} \left[ \int_0^T |Af(X^\varepsilon(s), Y^\varepsilon(s))|^q ds \right] < \infty.$

With $f \in C_b^2$, we can apply the Taylor-Lagrange formula :

$$ \sup_{t \le T} |R_f^\varepsilon(t)| \le \frac{\varepsilon}{2} \|f''\|_\infty \int_0^T X^\varepsilon(s) \left( R_0 + B_0 + D_0 \left[ 1 + (Y_{tot}^\varepsilon(s))^{m_0} \right] \right) ds, $$

given hypotheses 2.1, 2.2 et 2.4. Let $T^\varepsilon(t) = X^\varepsilon(t) + \sum_{i=1}^d Y_i^\varepsilon(t)$. From Proposition 3.1, there exists $C_1 > 0$ independent of $\varepsilon$ such that :

$$ \mathbb{E} \left( \sup_{t \le T} |R_f^\varepsilon(t)| \right) \le \frac{\varepsilon T}{2} \|f''\|_\infty C_1 \left( 1 + \sup_\varepsilon \sup_{s \le T} \mathbb{E} \left[ |T^\varepsilon(s)|^{m_0+1} \right] \right) \xrightarrow{\varepsilon \to 0} 0. $$

Also, from Proposition 3.1, for any $q > 1$, there exists, with similar calculations, $C_2 > 0$ independent of $\varepsilon$ such that :

$$ \mathbb{E} \left[ \int_0^T |Af(X^\varepsilon(s), Y^\varepsilon(s))|^q ds \right] \le \|f'\|_\infty^q C_2 T \left[ 1 + \sup_\varepsilon \sup_{s \le T} \mathbb{E} \left( |T^\varepsilon(s)|^{q(m_0+1)} \right) \right] \xrightarrow{\varepsilon \to 0} 0. $$

We turn to the compactness of $Y^\varepsilon$, which does not occur in a functional sense. Let $\Gamma_\varepsilon$ be the occupation measure of $Y^\varepsilon$ over time, defined by :

$$ \forall A \subset \mathbb{N}^d, \Gamma_\varepsilon([0, t] \times A) = \int_0^t 1_A(Y^\varepsilon(s)) ds. $$

It is clear that $\Gamma_\varepsilon \in \mathcal{L}_m(\mathbb{N}^d)$, the space of measure $\mu$ on $[0, \infty[ \times \mathbb{N}^d$ such that $\mu([0, t[ \times \mathbb{N}^d) = t$. For any $K > 0$, from Proposition 3.1 and the Markov inequality,

$$ \mathbb{P} \left( Y_{tot}^{\varepsilon}(t) \ge K \right) \le \frac{\sup_{\varepsilon} \sup_{t>0} \mathbb{E} (Y_{tot}^{\varepsilon}(t))}{K}. $$

For any $\nu > 0$, there exists a $K > 0$, such that, for all $t > 0$ and $\varepsilon > 0$,

$$ \mathbb{P} \left( Y_{tot}^{\varepsilon}(t) \le K \right) \ge 1 - \nu, $$

hence with $B_K = \{y \in \mathbb{N}^d, y_{tot} \le K\}$, $\inf_\varepsilon \mathbb{E}(\Gamma_\varepsilon([0, t] \times B_K)) \ge t(1-\nu)$, and we can conclude with [40, Lemma 1.3] that $(X^\varepsilon, \Gamma_\varepsilon)$ is relatively compact in $\mathcal{D}_{\mathbb{R}+} [0, \infty[ \times \mathcal{L}_m(\mathbb{N}^d)$.

## 4.2 Identification of the limit

Let $(x, \Gamma)$ be a limiting point (up to a sub-sequence) of $(X^\varepsilon, \Gamma_\varepsilon)$ in $\mathcal{D}_{\mathbb{R}+} [0, \infty[ \times \mathcal{L}_m(\mathbb{N}^d)$. As the first component $X^\varepsilon$ has vanishing jump sizes, namely (a.s.) for all $T > 0$,

$$ \sup_{t \in [0, T]} |X^\varepsilon(t) - X^\varepsilon(t^-)| \le \varepsilon \to 0, $$
$x$ is clearly a continuous function.

We start by identifying the limit measure $\Gamma$. Using definition 2.1 with $g$ depending only on the second variable, i.e. $\forall x \in \mathbb{R}^+$, $g(x, \cdot) = h(\cdot)$, we have

$$ 
\begin{aligned}
A^{\varepsilon}h(x, y) &= \frac{\lambda_0(y)}{\varepsilon}x[h(y + e_1) - h(y)] + \sum_{i=1}^{d-1} \frac{\lambda_i(y)}{\varepsilon}y_i[h(y + e_{i+1} - e_i) - h(y)] \\
&\quad + \sum_{i=1}^{d} \frac{\mu_i(y)}{\varepsilon}y_i[h(y - e_i) - h(y)] \\
&= \frac{1}{\varepsilon}L_xh(y),
\end{aligned}
$$

according to (1.2). By definition of the infinitesimal generator, the process

$$ \varepsilon[h(Y^{\varepsilon}(t)) - h(Y^{\varepsilon}(0))] - \int_{0}^{t} \int_{\mathbb{N}^{d}} L_{X^{\varepsilon}(s)}h(y)\Gamma_{\varepsilon}(ds, dy), \quad (1.15) $$

is a $\mathcal{F}_t^\varepsilon$-martingale. Applying [40, Lemma 1.5 d], since rates are polynomial (Hypothesis 2.4) and using Proposition 3.1, we have, for any $h \in C_b^2(\mathbb{N}^d)$,

$$ \lim_{\varepsilon \to 0} \int_0^t \int_{\mathbb{N}^d} L_{X^{\varepsilon}(s)}h(y)\Gamma_{\varepsilon}(ds, dy) = \int_0^t \int_{\mathbb{N}^d} L_{X(s)}h(y)\Gamma(ds, dy). $$

From the martingale (1.15), it follows that

$$ \int_0^t \int_{\mathbb{N}^d} L_{x(s)}h(y)\Gamma(ds, dy), $$

is a continuous martingale of bounded variations, hence it is identically zero, for any $h \in C_b^2(\mathbb{N}^d)$ and with probability one. Let $\nu$ be the measure defined by $\nu([0, t]) = \Gamma([0, t], \mathbb{N}^d)$ for all $t$. From [40, Lemma 1.4], there exists a process $\gamma$ such that, for any $h \in C_b^2(\mathbb{N}^d)$,

$$ \int_0^t \int_{\mathbb{N}^d} L_{x(s)}h(y)\Gamma(ds, dy) = \int_0^t \int_{\mathbb{N}^d} L_{x(s)}h(y)\gamma_s(dy)\nu(ds). $$

With the same arguments as in [40, Example 2.3], we have, almost surely,

$$ \int_{\mathbb{N}^d} L_{x(t)}h(y)\gamma_t(dy) = 0, $$

for all $t$ outside of a set of zero Lebesgue measure. Since $Y$ has a unique invariant probability measure (theorem 3.1), we get $\gamma_t = \pi_{x(t)}$, and $\Gamma(ds, dy) = \pi_{x(ds)}(dy) ds$.

It only remains to identify function $x$. With the martingale $M_f^\varepsilon(t)$ defined in (1.14), and with $f = Id$, taking the limit $\varepsilon \to 0$, there exists a filtration $(\mathcal{F}_t)_t$ such that (applying [40, Lemma 1.5 d] again) :

$$ M_t = x(t) - x^{in} - \int_0^t \int_{\mathbb{N}^d} AId(x(s), y)\pi_{x(ds)}(dy) ds, $$

is a $(\mathcal{F}_t)_t$-martingale. For $s > 0$,

$$ 
\begin{aligned}
\int_{\mathbb{N}^d} AId(x(s), y)\pi_{x(s)}(dy) &= x(s)\int_{\mathbb{N}^d} [r_0(y) - \lambda_0(y) - \mu_0(y)]\pi_{x(s)}(dy) \\
&= x(s)\Lambda_0(x(s)).
\end{aligned}
$$
Furthermore, according to (1.14), $M_{Id}^{\epsilon}$ is a $\mathcal{F}_{t}^{\epsilon}$-martingale with predictable quadratic variation

$$ \langle M_{Id}^{\epsilon} \rangle_t = \epsilon \int_0^t [r_0(Y^{\epsilon}(s)) + \lambda_0(Y^{\epsilon}(s)) + \mu_0(Y^{\epsilon}(s))] ds. $$

Using Proposition 3.1, $\mathbb{E}(\langle M_{Id}^{\epsilon} \rangle_t) \underset{\epsilon \to 0}{\longrightarrow} 0$. Hence, $M_t \equiv 0$ for all times $t$, and $x$ is solution of

$$ x(t) = x^{in} + \int_{0}^{t} \Lambda_{0}(x(s))x(s) \, ds. $$

Since $S_f$ is continuous for any $|f| \le c(V_p + 1)$ (see theorem 3.1), we deduce by composition that $x \in C^1(\mathbb{R}^+)$, and is solution of:

$$ \frac{dx}{dt}(t) = \Lambda_0(x(t))x(t), \quad x(0) = x^{in}. $$

## 4.3 Uniqueness of the limit

With theorem 3.1, we have proved that, for any limiting point $x$ of $X^{\epsilon}$, the limit of the fast process $Y^{\epsilon}$ is uniquely determined by the invariant probability measure $\pi_x$ of $L_x$. It remains to prove the uniqueness of $x$, as a solution of

$$ \frac{dz}{dt} = \Lambda_0(z)z, \quad z(0) = x^{in}. \qquad (1.16) $$

We will prove that the right-hand side function $F: \mathbb{R}^+ \to \mathbb{R}$ defined by

$$ F(x) = x\langle r_0 - \lambda_0 - \mu_0, \pi_x \rangle, $$

is locally Lipschitz continuous. Let $x_1 \in \mathbb{R}^+$ and $U \subset \mathbb{R}^+$ a neighborhood of $x_1$. From hypotheses 2.1 et 2.4, we get, for any $x_2 \in U$:

$$ \begin{aligned} |F(x_1) - F(x_2)| &= |\langle r_0 - \lambda_0 - \mu_0, \pi_{x_1} - \pi_{x_2} \rangle x_1 + \langle r_0 - \lambda_0 - \mu_0, \pi_{x_2} \rangle (x_1 - x_2)| \\ &\le x_1 |\langle r_0 - \lambda_0 - \mu_0, \pi_{x_2} - \pi_{x_1} \rangle| \\ &\quad + |x_1 - x_2|(R_0 + B_0 + D_0)\langle 1 + (y_{tot})^{m_0}, \pi_{x_2} \rangle. \end{aligned} $$

Choosing $p = (m_0+d-1, m_0+d-2, ..., m_0) \in \mathcal{P}$, $1+(y_{tot})^{m_0} \le 1+V_p$. From hypothesis 2.4 et theorem 3.1, there exists $C_1 > 0$ (depending on $x_1$) such that

$$ \forall x_2 \in U, \quad |F(x_1) - F(x_2)| \le C_1|x_1 - x_2|. $$

$F$ is thus locally Lipschitz continuous on $\mathbb{R}^+$, and we can deduce the uniqueness of a maximal solution of (1.16) on $\mathbb{R}^+$ from the Cauchy-Lipschitz theorem. Since $\Lambda_0$ is upper-bounded by $R_0$, we deduce from the comparison principle that $x$ is a global solution of (1.16).

# 5 Numerical illustration and biological application

In this section, we present some results on the linear version of the model (with constant rates), and on the identifiability of the model parameters. We then detail an algorithm suited for simulating $\pi_x$, which we use to illustrate the numerical convergence. Finally, we compare the model behavior in different parameter regimes specifically designed for our application to ovarian follicle development.
## 5.1 Linear case

Let us consider the process $(X^\epsilon(t), Y^\epsilon(t))_{t\ge0}$ (see definition 2.1) with all event rates $\lambda_i, \mu_i$ fixed to constant values. From theorem 2.1, we have:

**Proposition 5.1 (Linear case)**

For any initial condition $(X^\epsilon(0), Y^\epsilon(0)) \in (\epsilon\mathbb{N} \times \mathbb{N}^d)$ such that, for some $x^{in} \in \mathbb{R}^+$,

$$ \forall k > 1, \sup_{\epsilon} \mathbb{E}(X^{\epsilon}(0)^k) < \infty, \sup_{\epsilon} \mathbb{E}(Y_{tot}^{\epsilon}(0)^k) < \infty, X^{\epsilon}(0) \xrightarrow{a.s} x^{in} $$

the process $(X^\epsilon, Y^\epsilon)$ has a unique limit process $(x, \pi) \in \mathcal{C}^1(\mathbb{R}^+) \times \mathcal{L}_m(\mathbb{N}^d):$

$$ \begin{cases} x(t) = x^{in} \exp\left((r_0 - \lambda_0 - \mu_0)t\right), \\ \pi_{x(t)}(y) = \prod_{i=1}^{d} \frac{(x(t)\alpha_i)^{y_i}}{y_i!} e^{-x(t)\alpha_i} \end{cases} $$

with $\alpha_i = \prod_{j=0}^{i-1} \frac{\lambda_j}{\lambda_{j+1} + \mu_{j+1}}$ for $1 \le i \le d$.

**Proof:** The result follows directly from theorem 2.1, and from the forward Kolmogorov equation applied to the fast process Y [53, Thm 2.4].

Suppose now we have $\bar{x}(t)$ and $(\bar{y}_j(t))_{1\le j\le d}$ the number of individuals (e.g. follicles) observed in each compartment at any time between 0 and $T_f$. Since the number $x$ of individuals in the slow compartment (reserve) has an analytic expression, we can identify $\Lambda_0$ and $x^{in}$ from $(\bar{x}(t))_{0\le t\le T_f}$ for large $T_f$. Following [67], we use the following quantities to compute $x^{in}$ and $\Lambda_0$:

$$ \begin{aligned} \int_0^{T_f} \bar{x}(s) \, ds &= \frac{x^{in}}{\Lambda_0} (1 - \exp(-\Lambda_0 T_f)) \sim \frac{x^{in}}{\Lambda_0}, \\ \int_0^{T_f} s \bar{x}(s) \, ds &= \frac{x^{in}}{\Lambda_0^2} - \frac{x^{in}}{\Lambda_0} (T_f + \frac{1}{\Lambda_0}) \exp(-\Lambda_0 T_f) \sim \frac{x^{in}}{\Lambda_0^2}. \end{aligned} \quad (1.17) $$

To identify separately the exhaustion rate $\lambda_0 + \mu_0$ of the first compartment from the renewing rate $r_0$, we need additional observations to distinguish between outgoing and from births. To estimate the remaining parameters, $(\alpha_i)_{1\le i\le d}$, we apply the maximum likelihood principle to the Poisson distribution of each compartment $j$, for any $1 \le j \le d$ $(\bar{y}_j(t))_{0\le t\le T_f}$, and obtain an estimate $\bar{\alpha}_j$ of $\alpha_j$ as:

$$ \bar{\alpha}_j = \left( \int_0^{T_f} \bar{y}_j(t) \, dt \right) / \left( \int_0^{T_f} \bar{x}(t) \, dt \right). \quad (1.18) $$

To identify separately the $\lambda_i$s from the $\mu_i$s, additional observations would be needed on the number of deaths occurring in each compartment, or the transit times from one compartment to the next.
<figure><img src="image_13.png" /><figcaption>FIGURE 1.1 – Main steps of the limit model simulation</figcaption></figure>

## 5.2 Strategy to compute $\pi_x$ and the limit model

We use an Euler scheme over the time range $[0, 5]$ (with a time step of $5 \times 10^{-3}$) to simulate the ordinary differential equation ruling the dynamics of the first compartment. To evaluate the right hand side, the scheme requires the value of $\Lambda_0$ at each time step, given the current value of $x$ (Subsection 5.2).

Let $x \in \mathbb{R}^+$, $N_0 \in \mathbb{N}^*$ and $T_f > 0$. Algorithm 1 computes both the expectation $E$ of $\pi_x$, and the corresponding rate $\Lambda_0(x) = \sum_{y \in \mathbb{N}^d} (r_0(y) - \lambda_0(y) - \mu_0(y))\pi_x(y)$. The Markov chain associated with the fast process is first simulated $N_0$ times over the time range $[0, T_f]$ by running an exact simulation algorithm [28], then the computed trajectories (time and location $(t_i)_i, (Y_{t_i})_i$ at each jump in $[0, T_f]$) are used to estimate $\pi_x$. The values of $T_f$ and $N_0$ are a good compromise between accuracy (owing to [65], the error on $\lambda_0$ can be arbitrary small for large $T_f$ and $N_0$) and speed.

**Algorithm 1 Nested algorithm (x, T_f, N_0)**

$Y_{start} \leftarrow \text{Empty matrix of size } (N_0 + 1, d)$

$\Lambda_0 \leftarrow 0$

**for** n = 1,...,N_0 **do**

*(t_i)_i, (Y_{t_i})_i \leftarrow \text{Exact Simulation Algorithm starting in } Y_{start}[n] \text{ on } [0, T_f].*

$\Lambda_0 \leftarrow \Lambda_0 + \sum_i (r_0(Y_{t_i}) - \lambda_0(Y_{t_i}) - \mu_0(Y_{t_i}))(t_{i+1} - t_i)$

$E \leftarrow E + \sum_i Y_{t_i}(t_{i+1} - t_i)$

$Y_{start}[n+1] \leftarrow Y_{T_f}$

**end for**

**return** $\Lambda_0, E$

## 5.3 Parameter values

We specify the formulation of the event rates as proposed in [6] for a toy model of mammal follicle development. The numerical study performed here does not intend to reproduce accurately the follicle numbers in a specific species, and rather gives semi-quantitative insight into the follicle population dynamics. We focus on the control of the exit from the reserve (given $r_0 = 0$) by the growing follicles ($d = 3$), with $\lambda_0(y)$ expressed as :

$$ \lambda_0(y) = a + \frac{b}{1 + c(y_1 + y_2 + y_3)} \quad (1.19) $$
A schematic scheme of the 3d process $(X, Y) \in D_{\mathbb{R}^+ \times \mathbb{N}^3}[0, \infty[$ is provided in Figure 1.6. We instantiated the parameter values in $\lambda_0(y)$ so that the control exerted onto the reserve be operating (see Figure 1.7) : $a = 0.1$, $b = 1.5$, $c = 0.01$. All other rates are fixed to constant values $\lambda_1(y) = \lambda_2(y) = 0.5$, $\mu_0(y) = \mu_1(y) = \mu_2(y) = 0.3$, $\mu_3(y) = 0.2$, in such a way that the distribution $\pi_x$ of the limit model is non-monotonous (see Figure 1.8), as observed for growing follicles [63]. At initial time, the growing compartments are empty, and the reserve size is set to $x^{in} = 10^2$.

## 5.4 Illustration of the convergence

Subsection 5.4 illustrates the convergence of the sequence $(X^\epsilon, Y_1^\epsilon, Y_2^\epsilon, Y_3^\epsilon)$ to $(x, \pi)$ as $\epsilon \to 0$ (theorem 2.1). The trajectories of the initial, full model (definition 2.1) are obtained by averaging $10^4$ runs of an exact stochastic algorithm for different values of $\epsilon$, and compared to the limit model. To improve the estimate of the mean of $\pi_x$, we average the output of the Euler scheme combined with Algorithme 1 on $10^3$ runs. We quantify the convergence rate (top right panel) using the $l^1$ norm error at time $t$, $E_1$, and the total $l^1$ norm, $E_T$:

$$E_1(t) = |\mathbb{E}X_0^\epsilon(t) - x(t)| + \sum_{i=1}^3 |\mathbb{E}Y_i^\epsilon(t) - \mathbb{E}\langle x \mapsto x_i, \pi_t \rangle|, \quad E_T = \int_0^T E_1(t) \, dt.$$

<figure><img src="image_14.png" /><figcaption>FIGURE 1.2 – Convergence to the limit model</figcaption></figure>

Changes in the number $X^\epsilon$ of follicles in the reserve (top left panel) and in the number $(Y_i^\epsilon)_{1\le i\le 3}$ of growing follicles (bottom panels) for different values of $\epsilon$ (see legend insert). Top right panel : $l^1$ norm error $E_1(t)$ at time $t = 0.1$ (blue line) or $t = 2$ (orange line), and total $l^1$ norm $E_T$ (green line) as a function of $\epsilon$. The red dotted line provides a reference to an order 1 error in $\epsilon$.

The total $l^1$ norm error decreases linearly with $\epsilon$. This is also the case of the $l^1$ norm error after the transient regime (as for $E_1(2)$ at time $t = 2$), while the $\epsilon$-dependent pattern at earlier times (typically from time 0 to 0.5, as for $E_1(0.1)$ at time $t = 0.1$) is still decreasing, yet nonlinearly.

## 5.5 Biological model outputs

The status of the ovaries is mainly characterized by the number of quiescent follicles, as a function of age, on one side, and the distribution of growing follicles into different maturation stages, on the other side.
We now illustrate how the model outputs can participate in investigating two main issues : (i) What is the kinetics of the reserve decline up to (quasi-)exhaustion ?, and (ii) What does the growing follicle distribution look like and is its shape preserved with aging ?

We tackle the first question by studying the impact of the nonlinear formulation on the decay rate of the reserve, in comparison with linear formulations. To compare the linear and nonlinear formulations, we focus not only on the decrease in the number *x* of follicles in the reserve, but also on the cumulative number $E_3^T = \int_0^T \mu_3 Y_3(s) ds$ of follicles exiting from the last compartment. $E_3^T$ is a marker of performance of the reproduction process and can be assimilated to the number of ovulations occurring all along reproductive lifespan. To perform the comparison, we first simulate the nonlinear model with the same set of parameters as in subsection 5.2 to obtain synthetic data. We then fit the linear model to these data in three different ways (the parameter values are gathered in subsection 5.5). In the first two cases, we first estimate the meta-parameters $(\bar{\alpha}_i)_{1\le i\le 3}$ (1.18), and then the parameters affecting the reserve dynamics, $x^{in}$ and the decay rate $\lambda_0 + \mu_0$, using either a linear regression on $\ln(x)$, or the analytic expression of $x(t)$ (1.17). In the third case, we keep all parameters identical to the nonlinear formulation, except $\lambda_0$, whose value is retrieved from the least-squares fitting of the $(\bar{\alpha}_i)_{1\le i\le 3}$. Finally, we compute the expectation, variance and coefficient of variation of the normalized numbers of growing follicles ($(Y_i/x)_{i=1,2,3}$), numerically in the nonlinear case or analytically, using the intensities of the Poisson laws in Proposition 5.1 in the linear case).

<table><caption>**TABLE 1.2** – Parameter values used to compare the linear and non linear formulations of the model. Estimated values are in red. For the nonlinear model, the meta-parameters $(\bar{\alpha}_i)_{1\le i\le 3}$ are estimated using (1.18). These values are left unchanged for versions 1 and 2. For version 3, the value of $\lambda_0$ is retrieved from the least-squares fitting of the $(\bar{\alpha}_i)_{1\le i\le 3}$.</caption><thead><tr><th rowspan="3">Estimation method</th><th rowspan="3">Nonlinear model</th><th colspan="3">Linear models</th></tr><tr><th>version 1</th><th>version 2</th><th>version 3</th></tr><tr><th>linear regression on ln(x)</th><th>Integrals (1.17)</th><th>Estimation of λ<sub>0</sub> from (ᾱ<sub>i</sub>)<sub>1≤i≤3</sub></th></tr></thead><tbody><tr><td rowspan="3">(α<sub>i</sub>)<sub>1≤i≤3</sub></td><td>ᾱ<sub>1</sub> = 0.9490</td><td>ᾱ<sub>1</sub></td><td>ᾱ<sub>1</sub></td><td>α<sub>1</sub> = 0.9487</td></tr><tr><td>ᾱ<sub>2</sub> = 0.5931</td><td>ᾱ<sub>2</sub></td><td>ᾱ<sub>2</sub></td><td>α<sub>2</sub> = 0.5930</td></tr><tr><td>ᾱ<sub>3</sub> = 1.4821</td><td>ᾱ<sub>3</sub></td><td>ᾱ<sub>3</sub></td><td>α<sub>3</sub> = 1.4824</td></tr><tr><td>x<sup>in</sup></td><td>100</td><td>177.7</td><td>116.7</td><td>100.0</td></tr><tr><td>λ<sub>0</sub> + μ<sub>0</sub></td><td>0.1 + <sup>1.5</sup>&frasl;<sub>1+0.01(y<sub>1</sub>+y<sub>2</sub>+y<sub>3</sub>)</sub> + 0.3</td><td>1.532</td><td>1.237</td><td>0.758 + 0.3</td></tr></tbody></table>

The fitting results are illustrated on subsection 5.5. No linear version is able to capture the whole dynamics. Version 1 fits the reserve size well, yet there is a precocious shift on the cumulative number of follicles exiting from $Y_3$. In contrast, this output is rather well fitted by versions 1 and 2, yet at the expenses of a poorer fit to $x(t)$.

The effect of the nonlinear formulation is clearly seen on the pattern of the reserve dynamics. Indeed, the decay rate is not constant. It increases around the mid-simulation time, as a result of the decreasing number of follicles in the growing compartments, which progressively relaxes the feedback exerted onto the activation rate. Looking at the functional response given by the function $x \mapsto \Lambda_0(x)$, one has a direct access to the impact of growing follicles onto the reserve kinetics. In particular, the model results give mechanistic insight into statistical approaches that use nonlinear regression models to fit experimental data on the reserve size according to age in mammals (reviewed in [13], see also Figure 1.7).

The impact of nonlinear rates is further illustrated on Figure 1.4. The accelerated decay rate in
<figure><img src="image_15.png" /><figcaption>**FIGURE 1.3** – Comparison between the nonlinear model and three linear formulations fitted on simulated values of $x(t)$ and $Y_i(t)_{i=1,2,3}$.</figcaption></figure>

Each column represents the number of follicles in the reserve and fast compartments (from left to right). Each row illustrates the comparison of the nonlinear limit model (blue solid lines), with one of the three fitted linear models. In the first column, the vertical dashed lines represent the extinction time with a threshold level of one remaining follicle. See subsection 5.5 and text body for details on the fitting procedure and parameter values.

<figure><img src="image_16.png" /><figcaption>**FIGURE 1.4** – Expectation, variance and coefficient of variation of the normalized number of growing follicles as a function of $x$ in the linear and nonlinear cases.</figcaption></figure>

Blue lines: nonlinear formulation. Orange lines: linear formulations (version 3 has different values of $(\alpha_i)_{1\le i\le 3}$ (see subsection 5.5) that are not visually distinguishable from those of versions 1 and 2.). Each column corresponds to the normalized number of growing follicles in the successive compartments (from left to right: $Y_1/x$, $Y_2/x$, and $Y_3/x$).

the nonlinear case leads to an over-representation of growing individuals in the whole population for small reserve sizes.

In contrast, the distribution within the growing population $(Y_i/\sum_i Y_i)$ is conserved through time in both the nonlinear and linear cases (see Figure 1.8). This is a consequence of the convergence result on the shape of the growing follicle distribution. Indeed, for any time $t$, the distribution
has to be equal to the steady-state distribution associated with the generator, $L_{x(t)}$, which only evolves through the time-dependent input $x(t)$. The growing follicle distribution thus approximately preserves its shape through aging, while decreasing in amplitude, consistently with experimental observations in different species, [63]. Note also that the coefficient of variation is smaller than one in the nonlinear case, suggesting that the feedback mitigates the dispersion in individual numbers.

# 6 Conclusion

In this study, we have proved rigorously the convergence of a two-timescale stochastic model of population dynamics towards a reduced deterministic limit. The proof follows the general scheme of averaging techniques [40], that we had to adapt to our specific problem. We have obtained compactness and ergodicity properties thanks to a Foster-Lyapunov criterion [48] and coupling techniques. The main idea was to introduce a suitable change of variables to control the time evolution of the population. Due to the unidirectional motion from one compartment to the next, the population in each compartment can hardly be bounded separately, yet the cumulative sums of individual numbers can be controlled by a linear birth-migration process. This idea was used repeatedly to obtain proper control on the moments, the Foster-Lyapunov criterion, and irreducibility. We have also proved that the limit model has a unique solution, and thereby that the full sequence of rescaled processes converges to this solution. In more details, we had to prove that the quasi-stationary distribution of the fast component of the Markov chain is locally Lipschitz with respect to the slow component. To our knowledge, this issue is rarely tackled in details in works dedicated to stochastic averaging. To solve the uniqueness problem, we have introduced a suitable Poisson equation, a tool used in the field of sensitivity analysis. Finally, we have designed a numerical scheme dedicated to the limit model, and based on the nested algorithm [65], widely used in slow/fast problems in chemical kinetics.

From a biomathematical viewpoint, our study can be useful for compartmental systems or birth-death processes with unidirectional migration where the first compartment is much more crowded than the next ones, and has a slow dynamics. It is the case of developmental processes in which individuals in the most immature stage are both more abundant and subject to a slower dynamics than more mature individuals. Such processes are typically undergone in zoological groups spawning eggs and undergoing a metamorphosis process. Amongst these groups, some insect species have been the focus of population dynamics models (see for instance [60]), due to their role as disease vector or crop pest. The slow-fast hypothesis is natural in this context, and has even already been used to study a deterministic model of mosquito populations. In [61], the authors assume that “the egg stock is large, and its dynamics slow compared with the larvae stock”. In addition, they introduce a feedback of larvae onto the hatching term, so that, similarly as in our framework, the dynamics of the immature, slow stage is nonlinear and modulated by a faster stage.

From the application viewpoint, our nonlinear formulation generalizes previous attempts to model the dynamics of the ovarian follicle population along life. The limit model approach is based on a multiple timescale hypothesis. This hypothesis can be grounded on quantitative information to set the order of magnitude of the number of quiescent follicles, and of the activation rate. When monitoring the ovarian reserve, hence the number of ovarian follicles a woman (or female mammal) is endowed with at a given age, one usually separates the quiescent follicles from the growing ones. At birth, there are only quiescent follicles and the size of the reserve is approximately $5 \times 10^5$ (see [27] for a comparison between the human species and mouse). The total number of follicles continuously decreases until the time of menopause, tens of years later, when the follicle reserve is close to $10^3$
follicles. Moreover, while it takes in average fifty years to exhaust the pool of quiescent follicles, the whole growing sequence spans approximately three to four months, as established from cell kinetics studies [30] or grafts of ovarian cortex [17].

Linear population models were considered in a series of articles from MJ Faddy and co-authors, in mice and women (see e.g. [22] and [24]). Time-dependent coefficients were yet needed to achieve the parameter inference from follicle number data, revealing the underlying non-stationarity of the dynamics. Noteworthy, the parameter estimation from data-fitting leads to a much slower rate of activation of the first follicle class than the growth rate of the subsequent classes (corresponding to a small espilon value around $10^{-1}/10^{-2}$ according to the class number). The results provided in subsection 5.5 clearly illustrate that the nonlinear feedback terms impact the decay rate of the follicle reserve, which speeds up when the growing population diminishes, in consistency with the conclusions of statistical regression models (see e.g. the comprehensive studies [13, 64]). Many interesting and challenging questions remain to be investigated from a physiological and mathematical perspective. We intend to tackle first the issue of solving inverse problems fuelled by data on age-dependent follicle numbers in different species, and to reproduce experimental and physiological situations in which the reserve exhaustion occurs prematurely.
# 7 Appendix

## 7.1 Generator of the coupling

In the proof of lemma 3.1, we couple processes $(X^\epsilon, V^\epsilon)$ and $(M^\epsilon, W^\epsilon)$ with the following generator defined for all functions $h \in C_c^\infty(\mathbb{R}^+ \times \mathbb{N}^d)^2$:

$$ 
\begin{aligned}
C^\epsilon h(x, v, m, w) = & \frac{r_0(v)}{\epsilon} x\left(h(x + \epsilon, v, m + \epsilon, w) - h(x, v, m, w)\right) \\
& + \left(\frac{R_0}{\epsilon}m - \frac{r_0(v)}{\epsilon}x\right)^+ \left(h(x, v, m + \epsilon, w) - h(x, v, m, w)\right) \\
& + \frac{\lambda_0(v)}{\epsilon}x\left(h(x - \epsilon, v + \sum_{i=1}^d e_i, m, w + \sum_{i=1}^d e_i) - h(x, v, m, w)\right) \\
& + \left(\frac{B_0}{\epsilon}m - \frac{\lambda_0(v)}{\epsilon}x\right)^+ \left(h(x, v, m, w + \sum_{i=1}^d e_i) - h(x, v, m, w)\right) \\
& + \frac{1}{\epsilon}\left(\mu_0(v)x - \alpha_0 \min(x, m)\right)\left[h(x - \epsilon, v, m, w) - h(x, v, m, w)\right] \\
& + \frac{\alpha_0}{\epsilon}\min(x, m)\left[h(x - \epsilon, v, m - \epsilon, w + \sum_{i=1}^d e_i) - h(x, v, m, w)\right] \\
& + \frac{\alpha_0}{\epsilon}\left(m - \min(x, m)\right)\left[h(x, v, m - \epsilon, w + \sum_{i=1}^d e_i) - h(x, v, m, w)\right] \\
& + \sum_{i=1}^d \frac{\alpha_i u_i(v, w)}{\epsilon} \left[ \frac{\mu_i(v)}{\mu_i(v) + \lambda_i(v)}h\left(x, v - \sum_{j \ge i} e_j, m, w - e_i\right) + \frac{\lambda_i(v)}{\mu_i(v) + \lambda_i(v)}h\left(x, v - e_i, m, w - e_i\right) - h(x, v, m, w) \right] \\
& + \sum_{i=1}^d \frac{\alpha_i}{\epsilon}(w_i - w_{i-1} - u_i(v, w))\left[h\left(x, v, m, w - e_i\right) - h(x, v, m, w)\right] \\
& + \sum_{i=1}^d \frac{1}{\epsilon}\frac{\mu_i(v)}{\mu_i(v) + \lambda_i(v)}\left[(\mu_i(v) + \lambda_i(v))(v_i - v_{i-1}) - \alpha_i u_i(v, w)\right]\left[h\left(x, v - \sum_{j \ge i} e_j, m, w\right) - h(x, v, m, w)\right] \\
& + \sum_{i=1}^d \frac{1}{\epsilon}\frac{\lambda_i(v)}{\mu_i(v) + \lambda_i(v)}\left[(\mu_i(v) + \lambda_i(v))(v_i - v_{i-1}) - \alpha_i u_i(v, w)\right]\left[h\left(x, v - e_i, m, w\right) - h(x, v, m, w)\right]
\end{aligned}
$$
with $u_i(v, w) = \min(v_i - v_{i-1}, w_i - w_{i-1})$.

We verify that for any $x \le m$, and any function $h$ that depends only of $(x, v)$, we have $C^\epsilon h(x, v, m, w) = A^\epsilon h(x, v)$. Similarly, for any $h$ that depends only of $(m, w)$, $C^\epsilon h(x, v, m, w) = B^\epsilon h(m, w)$.

## 7.2 One additional lemma to upper bound $B^\epsilon V_p$ (lemma 3.2, subsection 3.1)

### Lemma 7.1

For all $i \in \mathbb{N}^*$ and $m \in \mathbb{N}$,
• $D_i^+(m) := (m+1)^i - m^i \le 2^i(m^{i-1} + 1),$

• $mD_i^-(m) := m[(m-1)^i - m^i] \le -m^i.$

**Proof:** For any $m \ge 1$:

$$
\begin{align*}
D_i^+(m) &= \sum_{j=0}^{i} \binom{i}{j} m^j - m^i = \sum_{j=0}^{i-1} \binom{i}{j} m^j \\
&= m^{i-1} \sum_{k=0}^{i-1} \binom{i}{k+1} \frac{1}{m^k} \\
&\le m^{i-1} \sum_{k=0}^{i-1} \binom{i}{k+1} \text{ since } m \ge 1 \\
&\le 2^i m^{i-1} \le 2^i (m^{i-1} + 1).
\end{align*}
$$

Given that $D_i^+(0) = 1 \le 2^i$, we get the first inequality.

$$
\begin{align}
m D_i^- (m) &= \sum_{j=0}^{i} \binom{i}{j} (-1)^j m^{i+1-j} - m^{i+1} = \sum_{j=1}^{i} \binom{i}{j} (-1)^j m^{i+1-j} \nonumber \\
&= m^i \sum_{j=1}^{i} \binom{i}{j} (-1)^j m^{1-j} \nonumber \\
&= m^i \underbrace{\sum_{k=0}^{i-1} \binom{i}{k+1} (-1)^{k+1} \frac{1}{m^k}}_{:=f_i(m)}. \tag{1.20}
\end{align}
$$

Since $f_i$ is continuous on $[1, \infty[$, $f_i(1) = -1$ and $\lim_{m \to \infty} f_i(m) = -i$, $f_i$ has a maximum $m_i$ on $[1, \infty[$ (which can be reached in infinity). More precisely, $m_i$ can be reached either in 1 or in $n_i \in [1, \infty[$ satisfying $f'_i(n_i) = 0$.

We prove by induction on $i \ge 1$ that $m_i$ is reached in 1 and that $f_i(1) = -1$.

• Initialization : $f_1(m) = -1 \le -1$ for any $m \ge 1$.

• Inductive step : Let $i \ge 1$. We assume that $f_i$ reaches its maximum $-1$ in 1.
The derivation of (1.20) gives :

$$
\frac{d D_{i+1}^{-}}{d m}(m) = m^{i-1}\left[ i f_{i+1}(m) + m f'_{i+1}(m) \right].
$$

By definition of $D_{i+1}^-$, we have

$$
(i + 1)D_i^-(m) = m^{i-1} \left[ i f_{i+1}(m) + m f'_{i+1}(m) \right].
$$
If $m_i$ is reached on $]1, \infty[$, $f'_{i+1}(m_{i+1}) = 0$, and we have

$$ 
\begin{aligned}
f_{i+1}(m_{i+1}) &= \frac{(i+1)D_i^-(m_{i+1})}{im_{i+1}^{i-1}} \\
&= \frac{(i+1)m_{i+1}^i}{im_{i+1}^{i-1}}f_i(m_{i+1}) \\
&\le \frac{(i+1)m_{i+1}}{i}f_i(m_i) \\
&\le -\frac{(i+1)m_{i+1}}{i} \\
&< -1
\end{aligned}
$$

Moreover, $f_{i+1}(1) = \sum_{k=0}^{i} \binom{i+1}{k+1} (-1)^{k+1} = -1$.

Hence, we have proved that for any $i \ge 1$, and for any $m \ge 1$,

$mD_i^-(m) \le -m^i$, and it is also true for $m = 0$. This gives the second inequality and achieves the proof of (1.20).

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

## 7.3 Additional Figures

<figure><img src="image_18.png" /><figcaption>FIGURE 1.5 – Path from $v = (1, 5, 5, 6)$ to $w = (3, 4, 6, 9)$</figcaption></figure>

Cumulative sums of the number of individuals in state $y$ (red solid line), $z$ (blue solid line) and $v^{(1)}$ (green solid line). The black oriented arrows represent $w_i - v_i$ for each compartment $i$. The first step consists in adding as many individuals as necessary to state $v$ (red lines), so that state $v^{(1)}$ (green lines) overcomes state $w$ (blue lines) in all compartments. The next steps consists in compensating for the numbers of individuals in each compartment, first through migration events and then through death events in the last compartment.
<figure><img src="image_19.png" /><figcaption>**FIGURE 1.6** – Scheme of the model with three compartments, ($X^\epsilon, Y_1^\epsilon, Y_2^\epsilon, Y_3^\epsilon$), used to perform the numerical illustrations in section 5.
All rates are fixed to constant values except $\lambda_0$ (see (1.19)).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**FIGURE 1.7** – $\Lambda_0$ (theorem 2.1) as function of $x$ computed with Algorithm 1.
Blue lines: nonlinear model simulated with Algorithm 1. Orange (resp. green, red) line: version 1 (resp. version 2, version 3) of the linear formulation. The light-blue dashed line is the graph of the function $x \to 0.491 + \frac{1.24}{1+0.0134x}$.</figcaption></figure>

Figure 1.7 represents the decay rate $\Lambda_0$ of the slow compartment in the nonlinear and linear scenarios. From (1.19), in the nonlinear case, this rate evolves between $a + \mu_0$ and $a + b + \mu_0$. Note that the decay rates inferred from the linear models are all within this range. In the nonlinear scenario, the functional response has a similar shape as the graph of the function $x \to a + \frac{b}{1+cx}$ (with $a, b, c$ positive real numbers), which has been proposed as a statistical regression model for the number of quiescent follicles according to age (see subsection 5.5 in the main text).
<figure><img src="image_21.png" /><figcaption>**FIGURE 1.8** – Distribution of the growing follicles at different times
Snapshots of the distribution at time $t = 0.5$ (left panels), $t = 1.5$ (middle panels) and time $t = 3$ (right panels).
Upper panels : numbers of growing follicles. Lower panels : normalized numbers within the growing population ($Y_i/\sum_i Y_i$). Blue bars : nonlinear model. Orange (resp. green, red) bars : version 1 (resp. version 2, version 3) of the linear formulation.</figcaption></figure>

Figure 1.8 illustrates the auto-similarity of the quasi-stationary distribution in the fast compart-
ments. Note the almost perfect preservation of shape in the distribution of normalized numbers,
and the decrease in the bar amplitudes, as time progresses, in the distribution of non-normalized
numbers.

**8 Simulation numérique de π**

Cette section vise à détailler les différentes méthodes envisagées pour simuler le modèle limite étudié précédemment ainsi qu’à expliquer la méthode choisie dans la sous-section 5.2.

Nous rappelons le modèle limite obtenu plus tôt dans ce chapitre (Théorème 2.1) :

$$
\left\{
\begin{array}{l}
x'(t) = \Lambda_0(x(t))x(t), \quad t > 0, \text{ avec } x(0) = x^{in} > 0 \\
\\
\Lambda_0(x(t)) = \sum_{y \in \mathbb{N}^d} (r_0(y) - \lambda_0(y) - \mu_0(y)) \pi_{x(t)}(y)
\end{array}
\right.
\qquad (1.21)
$$

où, $\pi_x$ est la mesure de probabilité stationnaire de la chaîne associée au générateur $L_x$ défini pour tout $x > 0$ et pour tout $f \in C_b(\mathbb{N}^d)$ par :

$$
\begin{equation}
\begin{split}
\forall y \in \mathbb{N}^d, L_x f(y) ={}& \lambda_0(y)x[f(y + e_1) - f(y)] + \sum_{i=1}^{d-1} \lambda_i(y)y_i[f(y + e_{i+1} - e_i) - f(y)] \\
& + \sum_{i=1}^{d} \mu_i(y)y_i[f(y - e_i) - f(y)].
\end{split}
\tag{1.22}
\end{equation}
$$

Ce modèle est composé de deux objets inter-dépendants de natures différentes : *x* la solution d'un système d'équations différentielles ordinaires (EDO) dépendant d'une moyennisation d'une fonction contre *π* ; et *π* la mesure invariante d'une chaîne de Markov dépendant de *x*. A priori
la simulation du modèle limite nécessite la mise en oeuvre simultanée d'une méthode pour chaque composante du modèle limite. Il existe de nombreux schémas numériques permettant de simuler une EDO (schéma d'Euler ou schémas de Runge-Kutta) et leurs propriétés (convergence, consistance et stabilité) sont bien connues. Pour cette raison la suite se focalisera sur la simulation de la mesure stationnaire $\pi$.

Plusieurs types de méthodes ont été développés afin de simuler la mesure invariante d'une chaîne de Markov [39]. La simulation du modèle limite consiste à chaque pas de temps $t_i$ à simuler $\pi_{x(t_{i-1})}$ afin de pouvoir résoudre numériquement l'équation différentielle ordinaire pour calculer $x(t_i)$. La méthode choisie pour simuler $\pi$ doit donc être peut couteuse en temps de calcul. Deux méthodes approchées sont présentées ci-dessous : une méthode de troncature de l'espace d'état (section 4.2 de [39]) et une méthode simulant la chaîne de Markov.

Afin d'alléger les notations, nous omettrons la dépendance en $x$ de la chaîne de Markov associée au générateur $L_x$ : soient $x$ fixé, $Y$ le processus de Markov de dimension $d$ càdlàg défini en loi par le générateur $L_x$, $Q$ sa matrice de transition et $\pi$ sa distribution stationnaire.

## 8.1 Troncature de l'espace d'état

L'espace d'état de la chaîne de Markov $Y$ est décomposé en une union disjointe d'ensembles finis (appelés "niveaux") :

$$ \mathbb{R}^d = \bigcup_{k=0}^{\infty} \mathcal{N}_k := \bigcup_{k=0}^{\infty} \{y \in \mathbb{R}^d, y^{\text{tot}} = k\}, $$

Comme les transitions possibles à partir d'un état du niveau $k$ permettent d'atteindre seulement le niveau précédent (mort dans l'un des compartiments), le niveau suivant (naissance dans le compartiment 1) ou de rester dans le même niveau (migration entre 2 compartiments); la matrice $Q$ est tridiagonale par blocs :

$$ Q = \begin{pmatrix} Q_0 & Q_0^0 & 0 & 0 & \dots \\ Q_1^- & Q_1^1 & Q_1^1 & 0 & \dots \\ 0 & Q_2^- & Q_2^2 & Q_2^2 & \dots \\ 0 & 0 & Q_3^- & Q_3^3 & \dots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix} $$

où la matrice $Q_+^k$ (resp. $Q_-^k$) décrit les transitions du niveau $k$ vers le niveau supérieur (resp. inférieur) et la matrice $Q_k^k$ les transitions au sein du niveau $k$.

$Y$ fait ainsi partie de la classe des level-dependent quasi-birth-death processes (LDQBDPs) pour lesquels des schémas d'approximation appropriés ont été étudiés (ref).

L'idée de la méthode de troncature est de se ramener à un espace d'état fini en approchant $\pi$ seulement sur un nombre fini de niveaux. La forme tridiagonale de $Q$ permet alors d'exprimer simplement cette approximation en fonction des différents blocs présents dans la matrice $Q$.

Ainsi pour une troncature donnée $S$, on cherche à approcher $\pi$ par la distribution conditionnelle $\pi(\cdot|S)$ :

$$ \pi(y|S) = \begin{cases} \frac{\pi(y)}{\pi(S)} & \forall y \in S \\ 0 & \forall y \notin S \end{cases}, \qquad (1.23) $$
Soit

$$S_l = \bigcup_{k=0}^{l-1} \mathcal{N}_k$$

la troncature obtenue en gardant seulement les l-premiers niveaux.

Comme $\mathcal{N}_0 = \{0\}$, pour tout $k \le l$, la restriction de la distribution conditionnelle au niveau $\mathcal{N}_k$ peut s'écrire :

$$\forall y \in \mathcal{N}_k, \pi(y|S_l) = \pi(0|S_l)T^k(0, y) \quad k = 0, 1, \dots, l-1 \quad (1.24)$$

avec pour tout $0 \le k \le l$, $T^k = \prod_{n=1}^k R^n$ où les matrices $(R^n)_{n \in \mathbb{N}}$ sont la solution minimale positive des équations suivantes (découlant de l'équation de Kolmogorov progressive)

$$\left\{ \begin{array}{l} Q_+^{n-1} + R^n Q^n + R^n R^{n+1} Q_-^{n+1} = 0 \quad \forall n \in \mathbb{N}^*. \\ Q^0 + R^1 Q_-^1 = 0 \end{array} \right. \quad (1.25)$$

La valeur en 0 de la distribution conditionnelle est définie afin de normaliser à 1 la masse totale de $\pi(\cdot|S_l)$ :

$$\pi(0|S_l) \times \sum_{k=0}^{l-1} \sum_{y \in \mathcal{N}_k} T^k(0, y) = 1$$

La difficulté principale est le calcul des matrices $(R^n)_{n \in [1,l]}$. Si l'on souhaite calculer $R^{n+1}$ connaissant $R^n$, le système (1.25) contient $|\mathcal{N}_{n-1}| \times |\mathcal{N}_n|$ équations pour $|\mathcal{N}_n| \times |\mathcal{N}_{n+1}|$ inconnues, ce qui rend la résolution impossible car la taille des niveaux $\mathcal{N}_k$ croît progressivement. Pour cette raison, les matrices $(R^n)_{n \in [1,l]}$ ne peuvent pas être calculés par récurrence ascendante. Par contre, sachant $R^{n+1}$, le système (1.25) a une unique solution :

$$R^n = -Q_+^{n-1} (Q^n + R^{n+1} Q_-^{n+1})^{-1} \quad (1.26)$$

Ainsi, les matrices $(R^n)_{n \in [1,l]}$ peuvent être calculées par récurrence descendante avec $(R^l)$ donnée. Différentes méthodes ont été proposées pour le choix de $(R^l)$ ([5, 7]). La plus simple consiste à approcher $R^l$ par la matrice nulle. Cette condition équivaut à empêcher toute transition du niveau $\mathcal{N}_l$ vers le niveau $\mathcal{N}_{l-1}$ ce qui est cohérent avec la troncature choisie pour approcher $\pi$. Une autre méthode possible consiste à approcher la matrice $R^l$ par une matrice $\tilde{R}^l$ en utilisant l'équation (1.26), puis à utiliser cette nouvelle matrice pour approcher $\pi$. Dans la suite, $R^l$ sera approchée simplement par la matrice nulle.

L'algorithme pour approcher $\pi$ par la distribution conditionnelle $\pi(\cdot|S_l)$ est détaillé ci-dessous en pseudo-code dans l'algorithme 2.

Cette méthode présente plusieurs avantages : elle nécessite peu de connaissances sur le processus Y, est simple à implémenter et a un faible coût de calcul. De plus, la distribution conditionnelle $\pi(\cdot|S_l)$ converge en variation totale vers la distribution stationnaire $\pi$ lorsque $l$ tend vers l'infini.

Lorsque l'espace d'état est de dimension quelconque, le nombre d'état présent au niveau $\mathcal{N}_k$ croît rapidement quand $k$ augmente. Le calcul et le stockage des matrices $(R^n)_{n \in [1,l]}$ deviennent alors très coûteux en mémoire et sont la principale limitation de cette méthode.
**Algorithm 2 Troncature de π**

$$ \pi \leftarrow (1, 0, 0, \dots, 0) $$

$$ R_l \leftarrow 0 $$

**for** i = 1 ... l-1 **do**

$$ R^{l-i} \leftarrow -Q_{+}^{l-i-1} (Q^{l-i} + R^{l-i+1} Q_{-}^{l-i+1})^{-1} $$

**end for**

**for** i = 1, ... l-1 **do**

$$ \pi(\mathcal{N}_i) \leftarrow \pi(\mathcal{N}_{i-1}) R^i $$

**end for**

**return** $\pi / \sum_y \pi(y)$

## 8.2 Méthode via la chaîne de Markov

Au lieu d'essayer d'approcher $\pi$ pour un ensemble d'états, cette méthode se base sur la structure du système (1.28). En effet, afin de pouvoir simuler l'équation différentielle régissant la dynamique de $x$, il suffit de pouvoir calculer la fonction $x \mapsto \Lambda_0(x)$. Pour ce faire, la méthode ci-dessous utilise la simulation de la chaîne de Markov ainsi que le théorème ergodique suivant :

**Théorème 8.1**

Soit $Z$ un processus de Markov à temps continu ayant pour espace d'état $\mathbb{R}^d$.

Supposons que $Z$ a une unique mesure stationnaire $\mu$, alors pour tout $f$ bornée,

$$ \frac{1}{T} \int_{0}^{T} f(Z_{s}) \, ds \underset{T \to \infty}{\longrightarrow} \int f \, d\mu $$

Soit $f = r_0 - \lambda_0 - \mu_0$ et $T_f > 0$. En simulant le processus $Y$ pendant un temps $T_f$ en partant de $y_0$, notons $(T_i)_{0 \le i \le N-1}$ les temps de saut du processus (avec $T_0 = 0$ et $T_N = T_f$). Alors d'après le théorème ergodique, $\Lambda_0$ peut être approché par

$$ \begin{aligned} \tilde{\Lambda}_0(x) &= \int f \, d\pi_x \\ &\approx \frac{1}{T_f} \int_0^{T_f} f(Y_s) \, ds \\ &\approx \frac{1}{T_f} \sum_{k=0}^{N-1} f(Y_{T_k}) (T_{k+1} - T_k) \end{aligned} $$

Afin de s'affranchir le plus possible de la condition initiale $y_0$, la procédure précédente est répé- tée $N_0$ fois en partant à chaque fois d'une condition initiale différente (souvent tirée aléatoirement). L'algorithme obtenu est le suivant :
**Algorithm 3 Nested algorithm (x, T<sub>f</sub>, N<sub>0</sub>)**

$Y_{start} \leftarrow \text{Matrice des conditions initiales de taille } (N_0 + 1, d)$

$\Lambda_0 \leftarrow 0$

for n = 1,...,N<sub>0</sub> do

$(t_i)_i, (Y_{t_i})_i \leftarrow \text{Algorithme de simulation exacte commençant en } Y_{start}[n] \text{ pendant un temps } T_f.$

$\Lambda_0 \leftarrow \Lambda_0 + \sum_i (r_0(Y_{t_i}) - \lambda_0(Y_{t_i}) - \mu_0(Y_{t_i}))(t_{i+1} - t_i)$

$E \leftarrow E + \sum_i Y_{t_i}(t_{i+1} - t_i)$

end for

return $\Lambda_0, E$

## 8.3 Illustration

Afin d'illustrer les 2 schéma numériques présentés ci-dessus, nous nous plaçons dans le cas d'un modèle avec une réserve X et 3 compartiments plus rapides ($Y_1, Y_2, Y_3$). Considérons les taux de mort et de migration suivants :

$$
\left\{
\begin{array}{l}
\lambda_0(y) = 0.1 + \frac{1.5}{1+0.05(y_1+y_2+y_3)} \\
\lambda_1(y) = \lambda_2(y) = 0.5 \\
\mu_i(y) = 0.3, \text{ pour } 1 \le i \le 3 \\
\mu_4(y) = 0.2
\end{array}
\right.
\qquad (1.27)
$$

Soit $x(0) = 10$, la Figure 1.9 illustre la dynamique du système jusqu'à la mort de tous les follicules où la distribution $\pi$ est calculée avec la méthode de troncature (Sous-section 8.1) et l'équation différentielle de la réserve avec un schéma d'Euler explicite.

Pour les compartiments rapides (panel en haut à droite et les 2 panels du milieu), il n'y a pas de représentation trajectorielle. L'effectif moyen de chaque loi marginale de la distribution $\pi$ est donc représenté à la place.

La Figure 1.9 montre aussi un exemple de la distribution $\pi$ au temps $t \sim 1$ lorsque $x(t) = 3$. L'espace d'états de la distribution de probabilité est ordonné selon l'ordre lexicographique et est divisé en niveaux (représentés par l'alternance des bandes grises et blanches).

La Figure 1.10 représente un zoom de la distribution $\pi_3$ sur 3 niveaux d'états. Au sein de chaque niveau, le même motif semble se répéter avec plus ou moins d'amplitude.

L'algorithme de troncature d'état permet d'obtenir des informations sur la forme de la distribution de probabilité. Par contre, il ne peut pas être utilisé pour les populations de follicules à cause de 2 limitations :

* lorsque le nombre de compartiments rapides augmente, le nombre d'état dans chaque niveau augment aussi et le calcul de $\pi$ nécessite de manipuler de très grandes matrices,

* lorsque la condition initiale de la réserve $x(0)$ augmente, l'algorithme ne termine plus en un temps raisonnable car chaque itération du schéma d'Euler nécessite le calcul de $\pi_x$.

L'algorithme nested n'étant pas impacté par ces limitations, nous l'avons utilisé pour illustrer le modèle limite dans la sous-section 5.5.
<figure><img src="image_22.png" /><figcaption>**FIGURE 1.9** – Dynamique du système calculé à l'aide de l'algorithme de troncature</figcaption></figure>

Panel en haut à gauche : évolution au cours du temps de la réserve. Panel en haut à droite et panels du milieu : effectifs moyens des lois marginales de $\pi_{x(t)}$.

Panel du bas : Distribution de probabilité $\pi_3$ calculé avec l'algorithme de troncature. Les états sont ordonnés selon l'ordre lexicographique et séparés par niveaux (bandes grises).

<figure><img src="image_23.png" /><figcaption>**FIGURE 1.10** – Zoom sur les niveaux $N_9$, $N_{10}$ et $N_{11}$ de la distribution $\pi$ de la Figure 1.9</figcaption></figure>
# 9 Exploitation biologique du modèle limite

Cette section vise à illustrer l'intérêt du modèle limite grâce à sa confrontation avec les données biologiques présentées dans le chapitre 3.

Nous rappelons le modèle limite obtenu plus tôt dans ce chapitre (Théorème 2.1) :

$$ egin{cases} x'(t) = \Lambda_0(x(t))x(t), & t > 0, \text{ avec } x(0) = x^{in} > 0 \\ \Lambda_0(x(t)) = \sum_{y \in \mathbb{N}^d} (r_0(y) - \lambda_0(y) - \mu_0(y))\pi_{x(t)}(y) \end{cases} \quad (1.28) $$

où, $\pi_x$ est la mesure de probabilité stationnaire de la chaîne associée au générateur $L_x$ défini pour tout $x > 0$ et pour tout $f \in C_b(\mathbb{N}^d)$ par :

$$ \forall y \in \mathbb{N}^d, L_x f(y) = \lambda_0(y)x[f(y + e_1) - f(y)] + \sum_{i=1}^{d-1} \lambda_i(y)y_i[f(y + e_{i+1} - e_i) - f(y)] \\ + \sum_{i=1}^{d} \mu_i(y)y_i[f(y - e_i) - f(y)]. \quad (1.29) $$

Les deux objets mathématiques distincts composant le modèle limite (la solution $x$ d'une équation différentielle ordinaire, couplée avec une distribution quasi-stationnaire $\pi_x$) peuvent être exploité d'un point de vue biologique :

* la lente exhaustion de la réserve est un facteur clé de la capacité de reproduction. En effet, suffisamment de follicules doivent s'activer afin de permettre la reproduction à un instant donné de la vie de l'individu, mais la réserve ne doit pas s'épuiser trop vite sous peine de perdre définitivement la capacité de se reproduire.

Le taux de décroissance de la réserve est donc central et est représenté dans le modèle limite par la fonction $\Lambda_0$.

* La distribution des follicules selon leur maturation à un instant donné de la vie est parfaitement représenté par la distribution $\pi_x$.

Afin d'illustrer l'exploitation biologique du modèle limite, nous utilisons les mêmes jeux de données que ceux présentés et utilisés en détail dans la sous-section 2.1 du chapitre 3 ([23, 25]). Ils fournissent le nombre de follicules ovariens (classés en 5 compartiments) durant l'ensemble de la vie (de la naissance jusqu'à 500 jours) chez la souris.

Afin d'être cohérent avec le modèle étudié dans le chapitre 3, nous considérons ainsi un modèle à 5 compartiments (une réserve et 4 compartiments de maturation) et nous nous focalisons sur la régulation du taux d'activation $\lambda_1$ de la réserve :

$$ \lambda_0(y) = m_0 + \frac{f_0}{1 + K_0(y_2 + y_3)^n} $$

avec $m_0 > 0$ et $f_0, K_0$ et $n$ des réels positifs.

Le modèle limite est simulé en couplant la méthode (présentée dans la sous-section 8.2) permettant de simuler la distribution $\pi$ avec un schéma d'Euler explicite pour calculer la solution $x$ de l'équation différentielle ordinaire.
Pour la réserve, la variance du modèle de régression est estimée, par contre pour les compartiments, la variance utilisée est celle la distribution $\pi$ :

$$ 
\begin{aligned}
S_0(t_j) &= x(t_j) + \mathcal{N}(0, \sigma_0^2) \\
S_i(t_j) &= \mathbb{E}_i(\pi_{x(t_j)}) + \mathcal{N}(0, \mathbb{V}_i(\pi_{x(t_j)}))
\end{aligned}
$$

où,

* $\mathbb{E}_i(\pi_{x(t_j)})$ est la composante numéro $i$ de l'espérance de $\pi_{x(t_j)}$,

* $\mathbb{V}_i(\pi_{x(t_j)})$ est la composante numéro $i$ de la variance de $\pi_{x(t_j)}$.

Le modèle limite décrit la dynamique stationnaire de la croissance folliculaire. Pour cette raison, l'estimation est effectuée à partir du pic de cellules dans les stades de maturation (soit entre 100 et 500 jours) avec la méthode des moindres carrés. Le résultat de l'estimation est présenté sur la **Figure 1.11**.

<figure><img src="image_24.png" /><figcaption>**FIGURE 1.11** – Nombres de follicules estimés dans chaque compartiment.
Les lignes oranges représentent le nombre de follicules dans chaque compartiment, calculés avec les paramètres obtenus par l'estimation. Cette estimation est prolongée sur l'intervalle $[0, 100]$ (lignes vertes) pour faciliter la comparaison avec les données.</figcaption></figure>

Les nombres de follicules estimés dans chaque compartiment sont cohérents avec les données expérimentales après 100 jours. Par contre, en complétant la courbe entre 0 et 100 jours, nous pouvons constater que l'estimation n'est pas bonne et ce principalement pour la réserve et le dernier compartiment.

Pour la réserve, comme la décroissance est significative principalement pendant les 100 premiers jours, les données disponibles après 100 jours ne permettent pas d'estimer correctement la
vitesse de décroissance de la réserve. Pour les compartiments suivants, la dynamique de la phase transitoire, c'est-à-dire l'augmentation du nombre de follicules au début de la vie, ne peut pas être décrite avec le modèle limite.

Pour ces raisons, le chapitre 3 est consacré à l'extension de notre modèle afin de pouvoir tenir compte de l'ensemble des données expérimentales ainsi que de deux jeux de données additionnelles.
<figure><img src="image_25.png" /></figure>
# Chapitre 2

## Modèle à trois échelles de temps

Dans ce deuxième chapitre, le modèle à deux échelles de temps du chapitre 1 est raffiné et nous prouvons sous des hypothèses plus fortes l'existence et l'unicité d'un modèle limite.

Nous caractérisons ensuite l'erreur entre le modèle et sa limite.

Travail réalisé avec Léa Popovic.

## Sommaire

<table><tr><td>1</td><td>Motivation</td><td>53</td></tr><tr><td>2</td><td>Modèle à trois échelles de temps</td><td>53</td></tr><tr><td>2.1</td><td>Changements d'échelle en taille et en temps</td><td>53</td></tr><tr><td>2.2</td><td>Notations</td><td>54</td></tr><tr><td>2.3</td><td>Modèle obtenu</td><td>55</td></tr><tr><td>2.4</td><td>Hypothèses</td><td>56</td></tr><tr><td>3</td><td>Résultats préliminaires</td><td>57</td></tr><tr><td>3.1</td><td>Contrôle uniforme de Ψ<sup>ε</sup></td><td>57</td></tr><tr><td>3.2</td><td>Comportement en temps long de l'échelle rapide</td><td>58</td></tr><tr><td>3.3</td><td>Étude de l'échelle intermédiaire</td><td>59</td></tr><tr><td>4</td><td>Modèle limite</td><td>62</td></tr><tr><td>4.1</td><td>Compacité</td><td>63</td></tr><tr><td>4.2</td><td>Identification de la limite</td><td>63</td></tr><tr><td>4.3</td><td>Unicité de la limite</td><td>67</td></tr><tr><td>5</td><td>Discussion du modèle limite</td><td>67</td></tr><tr><td>6</td><td>Étude de l'erreur d'approximation</td><td>68</td></tr></table>
<table><tr><td>6.1</td><td></td><td>68</td></tr><tr><td>6.1.1</td><td>Décomposition de l'erreur</td><td>69</td></tr><tr><td>6.1.2</td><td>Sélection de la fonction H<sup>e</sup></td><td>69</td></tr><tr><td>6.1.3</td><td>Variation quadratique</td><td>71</td></tr><tr><td>6.1.4</td><td>Passage à la limite</td><td>72</td></tr><tr><td>6.2</td><td>erreur d'approximation pour trois échelles de temps</td><td>73</td></tr><tr><td>6.2.1</td><td>Décomposition de l'erreur</td><td>74</td></tr><tr><td>6.2.2</td><td>Contrôle de l'erreur issue de l'échelle rapide</td><td>75</td></tr><tr><td>6.2.3</td><td>Contrôle de l'erreur issue de l'échelle intermédiaire</td><td>76</td></tr><tr><td>6.2.4</td><td>Variation quadratique</td><td>78</td></tr><tr><td>6.2.5</td><td>Passage à la limite</td><td>79</td></tr><tr><td>6.3</td><td>Discussion de la limite obtenue</td><td>80</td></tr></table>
# 1 Motivation

Dans le chapitre 1, nous avons introduit un modèle à 2 échelles de temps afin de rendre compte de la différence de vitesse entre la réserve de follicules quiescents et les follicules en croissance. Nous avons ainsi considéré que l'ensemble de la croissance des follicules (après activation) s'effectuait avec le même ordre de grandeur de vitesse. Biologiquement, cette hypothèse traduit correctement l'activation des follicules et le début du développement folliculaire mais n'est pas vérifiée pour les follicules en fin de croissance (phase pré-ovulatoire par exemple). En effet, chez la femme, la maturation complète d'un follicule dure à peu près 5 mois alors que la croissance précédant l'ovulation ne dure que quelques semaines. En plus de cette différence de vitesse de croissance, il y a aussi une différence d'abondance :

* après la puberté, la réserve contient environ $10^4$ follicules quiescents et se vide progressivement,

* le nombre de follicules dans les stades terminaux est de l'ordre de $10^2$,

* il y a seulement quelques follicules pré-ovulatoire pour finalement aboutir à un seul follicule au moment de l'ovulation.

Cette diminution de l'abondance et cette accélération de la vitesse de croissance ne sont pas des caractéristiques propres à la femme mais se retrouve de façon plus ou moins marquées chez tous les mammifères.

Le développement des follicules va donc de plus en plus vite et au fur et à mesure de moins en moins de follicules sont présents dans les stades de croissance. Dans ce chapitre, nous raffinons le modèle en ajoutant une troisième échelle de temps pour prendre en compte un peu plus finement les différences de vitesse parmi les follicules en croissance. Puis, en adaptant la méthode utilisée dans le chapitre 1, les différences d'échelle de temps nous permettent de prouver l'existence et l'unicité d'un modèle limite sous certaines hypothèses.

Nous présentons ensuite un travail toujours en cours qui consiste à caractériser l'erreur entre le modèle et sa limite. La méthode utilisée est d'abord appliqué au modèle à deux échelles de temps puis au modèle plus compliqué comportant trois échelles de temps. Enfin, nous terminons ce chapitre en discutant les hypothèses faites pour caractériser l'erreur ainsi qu'une possible utilisation pratique de celle-ci.

# 2 Modèle à trois échelles de temps

## 2.1 Changements d'échelle en taille et en temps

Les stades de croissance sont séparés en trois catégories selon leur vitesse : la réserve *X* de follicules quiescents, les premiers stades de croissance ($Y_i$)$_{1\le i\le d}$ et enfin les stades de croissance restants ($Z_j$)$_{1\le j\le l}$. Comme précédemment, en numérotant les compartiments de 0 (X) à *d* + *l* ($Z_l$), un follicule dans le compartiment *k* peut migrer vers le compartiment *k* + 1 (excepté pour le dernier compartiment $Z_l$) avec le taux $\tilde{\lambda}_k$ ou bien mourir au taux $\tilde{\mu}_k$. Un taux de renouvellement $\tilde{r}_0$ est ajouté pour la réserve afin de modéliser certaines espèces (par exemple les poissons) pour lesquelles la réserve est renouvelée au cours de la vie.

Les transitions du processus ainsi obtenu sont détaillées dans la Table 2.1.
<table><caption>TABLE 2.1 – Transitions du processus $(X(t), Y(t), Z(t))_{t≥0}$ in $\mathbb{N}^{1+d+l}$.</caption><thead><tr><th>Compartiment</th><th></th><th>Transition</th><th>Taux</th></tr></thead><tbody><tr><td rowspan="3">Réserve X</td><td>Naissance</td><td>(X, Y, Z) → (X + 1, Y, Z)</td><td>r̃<sub>0</sub>(X, Y, Z)X</td></tr><tr><td>Mort</td><td>(X, Y, Z) → (X − 1, Y, Z)</td><td>μ̃<sub>0</sub>(X, Y, Z)X</td></tr><tr><td>Activation</td><td>(X, Y, Z) → (X − 1, Y + e<sub>1</sub>, Z)</td><td>λ̃<sub>0</sub>(X, Y, Z)X</td></tr><tr><td rowspan="3">(Y<sub>i</sub>)<sub>1≤i≤d</sub></td><td>Mort<br>(i ∈ [1, d])</td><td>(X, Y, Z) → (X, Y − e<sub>i</sub>, Z)</td><td>μ̃<sub>i</sub>(X, Y, Z)Y<sub>i</sub></td></tr><tr><td>Maturation,<br>i ∈ [1, d − 1]</td><td>(X, Y, Z) → (X, Y − e<sub>i</sub> + e<sub>i+1</sub>, Z)</td><td>λ̃<sub>i</sub>(X, Y, Z)Y<sub>i</sub></td></tr><tr><td>Maturation</td><td>(X, Y, Z) → (X, Y − e<sub>d</sub>, Z + e<sub>1</sub>)</td><td>λ̃<sub>d</sub>(X, Y, Z)Y<sub>d</sub></td></tr><tr><td rowspan="2">(Z<sub>j</sub>)<sub>1≤j≤l</sub></td><td>Mort j ∈ [1, l]</td><td>(X, Y, Z) → (X, Y, Z − e<sub>j</sub>)</td><td>μ̃<sub>d+j</sub>(X, Y, Z)Z<sub>j</sub></td></tr><tr><td>Maturation,<br>j ∈ [1, l − 1]</td><td>(X, Y, Z) → (X, Y − e<sub>j</sub> + e<sub>j+1</sub>, Z)</td><td>λ̃<sub>d+j</sub>(X, Y, Z)Z<sub>j</sub></td></tr></tbody></table>

Avec un abus de notation, soit $e = (e_i)_{i∈[1,d]}$ la base canonique de $\mathbb{N}^d$ et $e = (e_i)_{i∈[1,l]}$ celle de $\mathbb{N}^l$. Pour simplifier l'écriture, on pose $\tilde{\lambda}_{d+l} \equiv 0$.

Pour prendre en compte les différences d'abondance entre les différents types de follicules, nous introduisons un paramètre $\varepsilon$ et la renormalisation suivante est appliquée au système $(X, Y, Z)$ :

* le nombre de follicule dans les stades terminaux $((Z_j)_{1≤j≤l})$ est choisi comme référence et aucune renormalisation n'est appliquée à ces stades,

* la réserve contient énormément de follicules, disons de l'ordre de $\frac{1}{\varepsilon}$ ainsi $X$ devient $\varepsilon X$,

* les stades intermédiaires ont donc une abondance entre $\frac{1}{\varepsilon}$ et 1, soit $Y$ devient $\varepsilon^p Y$ avec $0 < p < 1$.

Enfin, le temps $t$ est accéléré en étant remplacé par $\frac{t}{\varepsilon}$ et nous faisons l'hypothèse suivante :

### Hypothèse 2.1

Soient pour tout $1 ≤ i ≤ d$, $1 ≤ j ≤ l$ et pour tout $(x, y, z) ∈ \mathbb{R}^+ × (\mathbb{R}^+)^d × \mathbb{N}^l$,

* $\lambda_i(x, y, z) = \frac{1}{\varepsilon^p} \tilde{\lambda}_i(x/\varepsilon, y/\varepsilon^p, z)$

* $\mu_i(x, y, z) = \frac{1}{\varepsilon^p} \tilde{\mu}_i(x/\varepsilon, y/\varepsilon^p, z)$

* $\lambda_{i+d}(x, y, z) = \tilde{\lambda}_{i+d}(x/\varepsilon, y/\varepsilon^p, z)$

* $\mu_{i+d}(x, y, z) = \tilde{\mu}_{i+d}(x/\varepsilon, y/\varepsilon^p, z)$

* $r_0(x, y, z) = \frac{1}{\varepsilon} \tilde{r}_0(x/\varepsilon, y/\varepsilon^p, z)$

* $\lambda_0(x, y, z) = \frac{1}{\varepsilon} \tilde{\lambda}_0(x/\varepsilon, y/\varepsilon^p, z)$

* $\mu_0(x, y, z) = \frac{1}{\varepsilon} \tilde{\mu}_0(x/\varepsilon, y/\varepsilon^p, z)$

Les fonctions $r_0$, $(\lambda_i)_{0≤i≤d+l}$ et $(\mu_i)_{1≤i≤d+l}$ définies ci-dessus sont indépendantes de $\varepsilon$.

## 2.2 Notations

Dans la suite, nous utiliserons les mêmes notations que le chapitre 1 rappelées ci-dessous :

* $C_b(\mathbb{N}^d) = \{f : \mathbb{N}^d \to \mathbb{R} \text{ bornée}\}$,

* $\mathbb{F} = \mathbb{R}^+ × (\mathbb{R}^+)^d × \mathbb{N}^l,$
* $\mathcal{D}_{\mathbb{N}^d}[0, \infty[$ : l'ensemble des fonctions à valeurs dans $\mathbb{N}^d$ qui sont continues à droite et limitées à gauche (càdlàg) sur $[0, \infty[$,

* $\mathcal{L}_m(\mathbb{N}^d)$ : l'ensemble des mesures $\mu$ sur $[0, \infty[ \times \mathbb{N}^d$ tel que $\mu([0, t[ \times \mathbb{N}^d) = t$,

Les notations suivantes seront utilisées de manière répétée :

* Pour $Z$ un processus $d$-dimensionnel appartenant à $\mathcal{D}_{\mathbb{N}^d}[0, \infty[$, $Z_{tot} = \sum_{i=1}^d Z_i$,

* Pour $y$ un vecteur de $\mathbb{R}^d$ et pour $1 \le i \le d$, $y_{i \to d} = (y_i, y_{i+1}, \dots, y_d)$.

* $\mathcal{P} = \{p \in (\mathbb{N}^*)^l, \forall i \in [1, l-1], p_i > p_{i+1}\}$

* Pour tout $p \in \mathcal{P}$, rappelons la définition de la fonction $V_p \in \mathcal{C}^1(\mathbb{N}^l)$ :

$$V_p(z) = \sum_{i=1}^{l} \left( \sum_{k=1}^{i} z_k \right)^{p_i}$$

Afin d'alléger les notations des propositions et théorèmes ci-dessous, une séquence de triplet de variables aléatoires $(s^\epsilon = (x^\epsilon, y^\epsilon, z^\epsilon))_\epsilon$ appartient à l'ensemble $\mathbb{F}_0$ si et seulement s'il existe $x \in \mathbb{R}^+$ tel que :

$$\begin{array}{l}
\forall \epsilon > 0, s^{\epsilon} = (x^{\epsilon}, y^{\epsilon}, z^{\epsilon}) \in \epsilon \mathbb{N} \times \epsilon^p \mathbb{N}^d \times \mathbb{N}^l \\
\forall k > 1, \quad \left\{
\begin{array}{ll}
\displaystyle \sup_{\epsilon} \mathbb{E}((x^{\epsilon})^k) & < \infty, \\
\displaystyle \sup_{\epsilon} \mathbb{E}((y^{\epsilon}_{tot})^k) & < \infty, \\
\displaystyle \sup_{\epsilon} \mathbb{E}((z^{\epsilon}_{tot})^k) & < \infty. \\
\multicolumn{1}{c}{x^{\epsilon} \xrightarrow{a.s.} x.}
\end{array}
\right.
\end{array}$$

$\mathbb{F}_0$ est l'ensemble des conditions initiales qui seront admissibles pour obtenir un modèle limite avec une condition initiale correctement définie.

## 2.3 Modèle obtenu

Le processus $\Psi^\epsilon = (X^\epsilon, Y^\epsilon, Z^\epsilon)$ obtenu suite aux changements d'échelles de la sous-section 2.1 est défini de la façon suivante :
**Definition 2.1**

Le processus de Markov $(\Psi^\epsilon(t))_{t\ge0}$ càdlàg (continu à droite et limité à gauche) a le générateur suivant, pour toute fonction $g : \mathbb{F} \to \mathbb{R}$ bornée,

$$
\begin{align*}
A^\epsilon g(x, y, z) &= \frac{r_0(x, y, z)}{\epsilon} x[g(x + \epsilon, y, z) - g(x, y, z)] + \frac{\mu_0(x, y, z)}{\epsilon} x[g(x - \epsilon, y, z) - g(x, y, z)] \\
&\quad + \frac{\lambda_0(x, y, z)}{\epsilon} x[g(x - \epsilon, y + \epsilon^p e_1, z) - g(x, y, z)] + \sum_{i=1}^d \frac{\mu_i(x, y, z)}{\epsilon} y_i[g(x, y - \epsilon^p e_i, z) - g(x, y, z)] \\
&\quad + \sum_{i=1}^{d-1} \frac{\lambda_i(x, y, z)}{\epsilon} y_i[g(x, y + \epsilon^p (e_{i+1} - e_i), z) - g(x, y, z)] \\
&\quad + \frac{\lambda_d(x, y, z)}{\epsilon} y_d[g(x, y - \epsilon^p e_d, z + e_1) - g(x, y, z)] \\
&\quad + \sum_{i=1}^l \frac{\mu_{i+d}(x, y, z)}{\epsilon} z_i[g(x, y, z - e_i) - g(x, y, z)] \\
&\quad + \sum_{i=1}^{l-1} \frac{\lambda_{i+d}(x, y, z)}{\epsilon} z_i[g(x, y, z + e_{i+1} - e_i) - g(x, y, z)]
\end{align*}
$$

Soit $(\mathcal{F}_t^\epsilon)_{t \ge 0}$ la filtration canonique associée à $(\Psi^\epsilon(t))_{t \ge 0}$.

**2.4 Hypothèses**

Le modèle décrit précédemment est très général en particulier car aucune hypothèse n'a été faite sur la forme des taux de transition. Enfin de pouvoir caractériser le modèle limite, nous supposons que :

**Hypothèse 2.2 (Dépendance des taux)**

Chaque compartiment peut réguler seulement les compartiments précédents et lui-même, c'est-
à-dire pour tout $(x, y, z) \in \mathbb{F}$,

$$
\forall 1 \le i \le d, \left\{ \begin{array}{l} \lambda_i(x, y, z) = \lambda_i(y_{i \to d}, z) \\ \mu_i(x, y, z) = \mu_i(y_{i \to d}, z) \end{array} \right.
$$

$$
\forall 1 \le i \le l, \left\{ \begin{array}{l} \lambda_{d+i}(x, y, z) = \lambda_{d+i}(z_{i \to l}) \\ \mu_{d+i}(x, y, z) = \mu_{d+i}(z_{i \to l}) \end{array} \right. .
$$

De plus, afin de simplifier les relations entre les différentes échelles de temps, les taux des réactions agissant sur deux échelles différentes ainsi que $r_0$ et $\mu_0$ sont de la forme suivante :

$$
\left\{
\begin{array}{ll}
r_0(x, y, z) = r_0(y, z) & \mu_0(x, y, z) = \mu_0(y, z) \\
\lambda_0(x, y, z) = \lambda_0(y, z) & \lambda_d(x, y, z) = \lambda_d(z)
\end{array}
\right.
$$

Afin de s'assurer de l'existence du processus $\Psi^\epsilon$ en tout temps ainsi que de l'existence de tous les objets mathématiques intervenant dans le modèle limite, nous faisons l'hypothèse suivante sur
les taux des réactions liants deux échelles de temps différentes :

**Hypothèse 2.3 (Taux de renouvellement et de migration inter-échelles)**

$$
\begin{align*}
\exists R_0 > 0, \forall (y, z) \in (\mathbb{R}^+)^d \times \mathbb{N}^l, \quad & r_0(y, z) \le R_0, \\
\exists B_0 > 0, \forall (y, z) \in (\mathbb{R}^+)^d \times \mathbb{N}^l, \quad & \lambda_0(y, z) \le B_0, \\
\exists B_d > 0, \forall (y, z) \in (\mathbb{R}^+)^d \times \mathbb{N}^l, \quad & \lambda_d(y_d, z) \le B_d.
\end{align*}
$$

Afin d'empêcher l'accumulation d'individu au sein d'un compartiment, nous faisons une hypo-
thèse sur le taux de migration et de sortie de chaque compartiment :

**Hypothèse 2.4 (Taux de migration λ)**

$$
\forall i \in [0, d+l], \forall v \in \mathbb{F}, \lambda_i(v) > 0,
$$

**Hypothèse 2.5 (Taux de sortie)**

$$
\forall i \in [0, d+l], \alpha_i := \inf_{v \in \mathbb{F}} (\lambda_i(v) + \mu_i(v)) > 0.
$$

Enfin, nous faisons l'hypothèse technique suivante :

**Hypothèse 2.6 (Taux polynomiaux)**

Il existe $m_0 \in \mathbb{N}$ et $D_0 > 0$ tels que pour tout $(y, z) \in (\mathbb{R}^+)^d \times \mathbb{N}^l$,

$$
\begin{cases}
\forall i \in [1, d+l-1], & \lambda_i(y, z) \le D_0 [1 + (y_{tot})^{m_0} + (z_{tot})^{m_0}] \\
\forall i \in [0, d+l], & \mu_i(y, z) \le D_0 [1 + (y_{tot})^{m_0} + (z_{tot})^{m_0}]
\end{cases}
$$

On suppose de plus que le taux de sortie R de la réserve défini pour tout (y, z) ∈ (R⁺)ᵈ × N'l
par

$$
R(y, z) = \lambda_0(y, z) + \mu_0(y, z) - r_0(y, z)
$$

est une fonction localement lipschitzienne selon sa première variable y.

3 Résultats préliminaires

Dans cette section, nous présentons l'ensemble des résultats nécessaires pour démontrer l'exis-
tence et l'unicité du modèle limite dans la section 4.

**3.1 Contrôle uniforme de Ψ<sup>ε</sup>**

**Proposition 3.1**

Pour toute condition initiale $\Psi^\epsilon(0) \in \mathbb{F}_0$, nous avons :

$$
\forall k > 1, \forall T > 0, \quad \sup_{\epsilon} \mathbb{E} \left( \sup_{t \le T} |X^{\epsilon}(t)|^k \right) < \infty,
$$

$$
\forall k > 1, \quad \sup_{\epsilon} \sup_{t \ge 0} \mathbb{E} \left( |Y_{tot}^{\epsilon}(t)|^k \right) < \infty,
$$

$$
\forall k > 1, \quad \sup_{\epsilon} \sup_{t \ge 0} \mathbb{E} \left( |Z_{tot}^{\epsilon}(t)|^k \right) < \infty.
$$
Afin de prouver la **Proposition 3.1**, nous définissons un processus stochastique ($M^\epsilon, N^\epsilon, P^\epsilon$) dont les taux de transitions sont linéaires (**Définition 3.1**). Grâce à un couplage entre ce processus et $\Psi^\epsilon$, nous montrons que ($M^\epsilon, N^\epsilon, P^\epsilon$) est un processus majorant de $\Psi^\epsilon$.

**Definition 3.1**

Soit $(M^\epsilon, N^\epsilon, P^\epsilon)$ le processus de Markov défini en loi par le générateur $B^\epsilon$ suivant, pour toute fonction $g : \mathbb{F} \to \mathbb{R}$ bornée :

$$
\begin{align*}
B^\epsilon g(m, n, p) ={}& \frac{R_0}{\epsilon} m[g(m + \epsilon, n, p) - g(m, n, p)] + \frac{\alpha_0}{\epsilon} m[g(m - \epsilon, n + \epsilon^p e_1, p) - g(m, n, p)] \\
& + \frac{B_0}{\epsilon} m[g(m, n + \epsilon^p e_1, p) - g(m, n, p)] \\
& + \sum_{i=1}^{d-1} \frac{\alpha_i}{\epsilon} n_i [g(m, n + \epsilon^p (e_{i+1} - e_i), p) - g(m, n, p)] + \frac{\alpha_d}{\epsilon} n_d [g(m, n - \epsilon^p e_d, p + e_1) - g(m, n, p)] \\
& + \frac{B_d}{\epsilon} n_d [g(m, n, p + e_1) - g(m, n, p)] \\
& + \sum_{i=1}^{l} \frac{\alpha_{i+d}}{\epsilon} p_i [g(m, n, p + e_{i+1} - e_i) - g(m, n, p)]
\end{align*}
$$

La preuve de la **Proposition 3.1** ne sera pas détaillée ici car elle utilise le même schéma de
preuve que la **sous-section 3.1 du chapitre 1** :

* changement de variables pour exprimer les composantes $y$ et $z$ en fonction des sommes cumulées $\left(\sum_{k=1}^{i} y_k\right)$ et $\left(\sum_{k=1}^{j} z_k\right)$, $1 \le i \le d$ et $1 \le j \le l$,

• contrôle des moments du processus linéaire ($M^\epsilon, N^\epsilon, P^\epsilon$),

• propagation du contrôle des moments grâce au couplage des processus $\Psi^\epsilon$ et ($M^\epsilon, N^\epsilon, P^\epsilon$).

**3.2 Comportement en temps long de l'échelle rapide**

Pour une valeur fixé de $Y^\epsilon$, nous introduisons et étudions le processus $Z$ représentant l'échelle de temps la plus rapide de $\Psi^\epsilon$. Pour la suite, soit $y \in (\mathbb{R}_*^+)^d$ un vecteur fixé.

**Proposition 3.2**

Soit $Z$ le processus de Markov càdlàg défini en loi, pour toute fonction $h \in C_b(\mathbb{N}^l)$, par le générateur $L_{y_d}$ suivant

$$
\begin{equation}
\forall z \in \mathbb{N}^l, L_{y_d} f(z) = \lambda_d(z) y_d[f(z + e_1) - f(z)] + \sum_{i=1}^{l-1} \lambda_{d+i}(z) z_i[f(z + e_{i+1} - e_i) - f(z)] \\
\phantom{\forall z \in \mathbb{N}^l, L_{y_d} f(z) = \lambda_d(z) y_d[f(z + e_1) - f(z)]} + \sum_{i=1}^{l} \mu_{d+i}(z) z_i[f(z - e_i) - f(z)]. \tag{2.1}
\end{equation}
$$

Alors $Z$ est exponentiellement ergodique et a une unique mesure de probabilité invariante $\pi_{y_d}$.
De plus, pour toute fonction $|f| \le c(1 + V_p)$ avec $c > 0$ et $p \in \mathbb{P}$, la fonction $S_f$ définie par
$$ \forall w \in \mathbb{R}^+, S_f(w) = \langle f, \pi_w \rangle, $$

est localement lipschitzienne sur $\mathbb{R}^+$.

**Démonstration:** Le processus $Z$ est similaire au processus $Y$ défini et étudié dans la *Sous-section 3.2* du chapitre 1. En utilisant le même schéma de preuve (preuve de l'irréductibilité de $Z$ puis existence d'une fonction de Lyapunov), on obtient l'unicité de $\pi_{y_d}$. L'introduction de l'équation de Poisson permet de conclure sur les propriétés de la fonction $S_f$.

## 3.3 Étude de l'échelle intermédiaire

Soient, pour tout $1 \le i \le d$ et pour tout $y \in (\mathbb{R}^+)^d$ :

$$ 
\left\{ 
\begin{aligned}
\Lambda_0(y) &= \int_{\mathbb{N}^l} \lambda_0(y, z) \pi_{y_d}(dz), & \Lambda_i(y) &= \int_{\mathbb{N}^l} \lambda_i(y, z) \pi_{y_d}(dz) \\
M_0(y) &= \int_{\mathbb{N}^l} \mu_0(y, z) \pi_{y_d}(dz), & M_i(y) &= \int_{\mathbb{N}^l} \mu_i(y, z) \pi_{y_d}(dz) \\
R_0(y) &= \int_{\mathbb{N}^l} r_0(y, z) \pi_{y_d}(dz)
\end{aligned}
\right. 
$$

Avec l'hypothèse faite sur les dépendances des taux de transition (Hypothèse 2.2), tous les taux définis ci-dessus dépendent seulement de la partie intermédiaire $y$.

Les compartiments les plus rapides interviennent sur la dynamique des compartiments intermédiaires via une moyennisation qui est matérialisée par les taux de transition ci-dessus. Afin de caractériser la limite du processus $\Psi^\epsilon$, nous avons besoin d'étudier le système d'EDOs associé à la partie intermédiaire du processus. Pour la suite, soit $x > 0$ un réel strictement positif.

**Definition 3.2**

Soit le système d'EDOs suivant :

$$ 
\left\{ 
\begin{aligned}
v'_1(t) &= \Lambda_0(v(t))x - (M_1(v(t)) + \Lambda_1(v(t)))v_1(t) \\
v'_i(t) &= \Lambda_{i-1}(v(t))v_{i-1}(t) - (M_i(v(t)) + \Lambda_i(v(t)))v_i(t), \text{ pour } 2 \le i \le d
\end{aligned}
\right. 
\quad (2.2) 
$$

et l'opérateur différentiel associé $\tilde{\mathcal{H}}_x$ défini pour toute fonction $g \in C^1(\mathbb{R}^d)$,

$$ \tilde{\mathcal{H}}_x g(y) = \Lambda_0(y) x \frac{\partial g}{\partial e_1}(y) - \sum_{i=1}^{d} (M_i(y) + \Lambda_i(y)) y_i \frac{\partial g}{\partial e_i}(y) + \sum_{i=1}^{d-1} \Lambda_i(y) y_i \frac{\partial g}{\partial e_{i+1}}(y) $$

Pour s'assurer que le système d'EDOs (2.2) a un unique point fixe, nous faisons l'hypothèse suivante sur la monotonie des taux de transition :
**Hypothèse 3.1 (Monotonie des taux)**

$\forall 0 \le i \le d-1,$ &nbsp;&nbsp;&nbsp;&nbsp; $\Lambda_i$ est décroissante selon $y_{i+1}, \dots, y_d,$
$\forall 1 \le i \le d-1,$ &nbsp;&nbsp;&nbsp;&nbsp; $\Lambda_i y_i$ est croissante selon $y_i,$
$\forall 1 \le i \le d,$ &nbsp;&nbsp;&nbsp;&nbsp; $(\Lambda_i + M_i)y_i$ est croissante selon $y_i, \dots, y_d.$

**Proposition 3.3**

Avec l'hypothèse 3.1,

le système (2.2) a un unique point fixe $v_x^* \in (\mathbb{R}^+)^d$.

**Démonstration :** Pour $1 \le i \le d$, soit $D_i = \Lambda_i + M_i$. Le système (2.2) peut s'écrire :

$$ egin{cases} v_1'(t) = F_1(v(t)) \\ v_i'(t) = F_i(v_{i-1 \to d}(t)) \text{ pour tout } 2 \le i \le d \end{cases} \qquad (2.3) $$

$$ \text{avec } \begin{cases} F_1(v) = \Lambda_0(v)x - v_1 D_1(v) \\ F_i(v_{i-1 \to d}) = \Lambda_{i-1}(v_{i-1 \to d})v_{i-1} - v_i D_i(v_{i \to d}) \text{ pour tout } 2 \le i \le d \end{cases}. $$

Nous allons construire un point fixe par récurrence, c'est-à-dire une solution $y^* \in (\mathbb{R}^+)^d$ de

$$ \begin{cases} F_1(y) = 0 \\ F_i(y_{i-1 \to d}) = 0 \text{ pour tout } 2 \le i \le d \end{cases} \qquad (2.4) $$

Prouvons par récurrence descendante pour $1 \le i \le d-1$ que

$$ \exists! \phi_i \text{ croissante telle que, } \forall y_d > 0, F_{i+1}(\phi_i(y_d), \dots, \phi_{d-1}(y_d), y_d) = 0 \qquad (2.5) $$

* **Initialisation** : Comme la dernière équation sur $F_d$ dépend seulement de $y_{d-1}$ et de $y_d$, et avec l'hypothèse sur la monotonie des taux (Hypothèse 3.1),

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial F_d}{\partial y_{d-1}} = \frac{\partial (\Lambda_{d-1} y_{d-1})}{\partial y_{d-1}} \ge 0 \\ F_d(0, y_d) < 0 \\ F_d(y_{d-1}, y_d) \ge \alpha_{d-1} y_{d-1} - y_d D_d(y_d) \underset{y_{d-1} \to \infty}{\longrightarrow} \infty \end{array} \right. $$

En appliquant le théorème des valeurs intermédiaires, nous avons le résultat voulu

$$ \forall y_d \ge 0, \exists! z = \phi_{d-1}(y_d) \ge 0 \text{ tel que } F_d(z, y_d) = 0 $$

Soient $y_d > 0$ et $h > 0$, avec l'hypothèse 3.1,

$$ \begin{align*} F_d(\phi_{d-1}(y_d), y_d + h) &= \Lambda_{d-1}(\phi_{d-1}(y_d), y_d + h)\phi_{d-1}(y_d) - (y_d + h)D_d(y_d + h) \\ &\le \Lambda_{d-1}(\phi_{d-1}(y_d), y_d)\phi_{d-1}(y_d) - y_d D_d(y_d) \\ &\le 0 \end{align*} $$

Comme $F_d$ est croissante selon $y_{d-1}$, $\phi_{d-1}(y_d + h) \ge \phi_{d-1}(y_d)$. En effectuant le même raisonnement pour $h$ négatif, la fonction $\phi_{d-1}$ est croissante.
• **Hérédité** : Soit $1 \le i \le d-2$ et supposons la propriété de récurrence (2.5) vraie pour tout $i+1 \le j \le d-1$.

Avec $\phi(y_{i+1 \to d}) = (\phi_{i+1}(y_d), \dots, \phi_d(y_d))$ ($\phi_d = I_d$),

$$F_{i+1}(y_i, \phi(y_{i+1 \to d})) = \Lambda_i(y_i, \phi(y_{i+1 \to d}))y_i - \phi_{i+1}(y_d)D_{i+1}(\phi(y_{i+1 \to d}))$$

Comme précédemment, avec l'hypothèse 3.1,

$$\left\{ \begin{array}{l} \frac{\partial F_{i+1}}{\partial y_i} = \frac{\partial}{\partial y_i} [y_i \Lambda_i(y_i, \phi(y_{i+1 \to d}))] \ge 0 \\ F_{i+1}(0, \phi(y_{i+1 \to d})) < 0 \\ F_{i+1}(y_i, \phi(y_{i+1 \to d})) \ge \alpha_i y_i - \phi_{i+1}(y_d) D_{i+1}(\phi(y_{i+1 \to d})) \underset{y_i \to \infty}{\longrightarrow} \infty \end{array} \right.$$

et en appliquant le théorème des valeurs intermédiaires

$$\forall y_d \ge 0, \exists! z = \phi_i(y_d) \ge 0 \text{ tel que } F_{i+1}(z, \phi_{i+1}(y_d), \dots, \phi_d(y_d)) = 0$$

Soient $y_d > 0, h > 0$. Pour simplifier le calcul ci-dessous, soient $z = (\phi_{i+1}(y_d), \dots, \phi_d(y_d))$ et $z_h = (\phi_{i+1}(y_d + h), \dots, \phi_d(y_d + h))$.

Comme les fonctions $\phi_k$ sont croissantes pour tout $k \ge i + 1$ et avec l'hypothèse 3.1,

$$\begin{align*} F_{i+1}(\phi_i(y_d), z_h) &= \Lambda_i(\phi_i(y_d), z_h)\phi_i(y_d) - \phi_{i+1}(y_d + h)D_{i+1}(z_h) \\ &\le \Lambda_i(\phi_i(y_d), z)\phi_i(y_d) - \phi_{i+1}(y_d)D_{i+1}(z) \\ &\le 0 \end{align*}$$

Comme $F_{i+1}$ est croissante selon $y_i$, $\phi_i(y_d + h) \ge \phi_i(y_d)$, c'est-à-dire que la fonction $\phi_i$ est croissante.

Il reste à identifier les valeurs possibles pour $y_d$ avec la seule équation que nous n'avons pas encore exploitée : $F_1 = 0$.

Soit $g$ la fonction définie sur $\mathbb{R}_*^+$ par $g(z) = F_1(\phi_1(z), \dots, \phi_{d-1}(z), z)$. Alors, pour tout $z > 0$,

$$\begin{align*} g' &= x \sum_{k=1}^{d} \frac{\partial \Lambda_0}{\partial y_k} \phi_k' - \phi_1' D_1 - \phi_1 \sum_{k=1}^{d} \phi_k' \frac{\partial D_1}{\partial y_k} \\ &\le 0 \end{align*}$$

Comme $\phi_1(0) = \dots = \phi_{d-1}(0) = 0$,

$$g(0) = \Lambda_0(0)x - 0 > 0$$

Avec l'hypothèse 3.1,

$$g(z) \le B_0 x - \phi_1(z) D_1(z) \underset{\|z\| \to \infty}{\longrightarrow} -\infty$$

D'après le théorème des valeurs intermédiaires, il existe un unique $z > 0$ tel que $g(z) = 0$. Par construction, $y^* = (\phi_1(z), \dots, \phi_{d-1}(z), z)$ est l'unique solution dans $(\mathbb{R}^+)^d$ du système (2.4).

Avec les hypothèses faites jusqu'ici, nous avons démontré l'existence et l'unicité d'un point fixe pour le système d'EDOs (2.2). Cependant pour garantir l'unicité de la limite du processus $\Psi^\epsilon$, il est nécessaire d'exclure les cas pour lesquels il y a existence d'un ou de plusieurs cycles limites. Pour des raisons techniques lors de la preuve (voir discussion en section 5), nous faisons l'hypothèse suivante sur les mesures stationnaires de l'opérateur $\tilde{\mathcal{H}}_x$ :
**Hypothèse 3.2**

Pour tout $x > 0$, il existe une unique mesure stationnaire de probabilité $\nu_x$ telle que, pour toute fonction $h : \mathbb{R}^d \to \mathbb{R}$ :

$$ \int_{\mathbb{R}^d} \tilde{\mathcal{H}}_x h(y) \nu_x(dy) = 0 \qquad (2.6) $$

Comme la mesure $\nu_x = \delta_{v_x^*}(dy)$ est solution de l'équation (2.6), avec l'hypothèse ci-dessus, c'est la seule mesure solution.

Pour garantir l'unicité du modèle limite, nous faisons aussi l'hypothèse suivante :

**Hypothèse 3.3**

Soit $\varepsilon > 0$, pour tout $x \in \mathbb{R}^+$, il existe $C_1 > 0$ (pouvant dépendre de $x$) tel que

$$ |v_{x+\varepsilon}^* - v_x^*| \le \varepsilon C_1 $$

Cette hypothèse semble raisonnable étant donné que l'opérateur $\tilde{\mathcal{H}}_x$ est linéaire en $x$ (voir la discussion en Section 5).

# 4 Modèle limite

Soit $\Gamma_\varepsilon$ la mesure d'occupation du couple $(Y^\varepsilon, Z^\varepsilon)$ au cours du temps :

$$ \forall A \subset (\mathbb{R}^+)^d, \forall B \subset \mathbb{N}^l, \Gamma_\varepsilon([0, t], A, B) = \int_0^t \mathbb{1}_{Y^\varepsilon(s) \in A, Z^\varepsilon(s) \in B} ds $$

Nous pouvons alors énoncer l'un des théorèmes principaux de ce chapitre :

**Théorème 4.1**

Pour toute condition initiale $\Psi^\varepsilon(0) \in \mathbb{F}_0$ (avec $x^{in}$ la limite de $X^\varepsilon(0)$), le processus $(X^\varepsilon, \Gamma_\varepsilon)$ est relativement compact dans $\mathcal{D}_{\mathbb{R}^+}[0, \infty[ \times \mathcal{L}_m(\mathbb{R}^d \times \mathbb{N}^l)$, et a un unique processus limite $(x, \Gamma) \in \mathcal{C}^1(\mathbb{R}^+) \times \mathcal{L}_m(\mathbb{R}^d \times \mathbb{N}^l)$ tel que :

$$ \left\{ \begin{array}{l} x'(t) = D_0(x(t))x(t), \ t > 0, \ \text{with } x(0) = x^{in} \\ D_0(x) = \sum_{z \in \mathbb{N}^l} (r_0(y_x^*, z) - \lambda_0(y_x^*, z) - \mu_0(y_x^*, z)) \pi_{(y_x^*)_d}(z) \\ \Gamma(dy, dz, ds) = \delta_{y_{x(s)}^*}(dy) \pi_{(y_{x(s)}^*)_d}(dz) ds \end{array} \right. $$

où,

* pour tout $x > 0$, $y_x^*$ est l'unique point fixe défini par la Proposition 3.3,

* pour $y \in \mathbb{R}^d$, $\pi_{y_d}$ est l'unique mesure de probabilité invariante définie par la Proposition 3.2.

Le preuve du Théorème 4.1 suit les mêmes étapes que la preuve dans le cas de 2 échelles de temps (Théorème 2.1 dans le chapitre 1) : un argument de compacité pour l'existence de la limite, puis l'identification de la limite et enfin l'unicité de cette limite.
La compacité du couple $(X^\epsilon, \Gamma_\epsilon)$ est obtenue grâce aux contrôles sur les moments (Proposition 3.1). L'identification de la limite nécessite de décomposer la mesure d'occupation limite afin d'identifier séparément la limite de l'échelle intermédiaire et celle de l'échelle rapide. Ce dernier point est la principale différence avec la preuve dans le cas de 2 échelles de temps.

## 4.1 Compacité

Commençons par la relative compacité de $X^\epsilon$ dans $\mathcal{D}_{\mathbb{R}^+}[0, \infty[$.

Nous avons fait les mêmes hypothèses que dans le chapitre 1 et ainsi obtenu le même contrôle sur les moments de $X^\epsilon$. Avec le même argumentaire, $X^\epsilon$ est relativement compact dans $\mathcal{D}_{\mathbb{R}^+}[0, \infty[$.

Passons maintenant à la compacité de la mesure d'occupation $\Gamma_\epsilon$ dans $\mathcal{L}_m(\mathbb{R}^d \times \mathbb{N}^l)$.
D'après l'inégalité de Markov et le contrôle sur les moments (Proposition 3.1), pour tout $K > 0$,

$$ \mathbb{P}(Y_{tot}^{\epsilon}(t) + Z_{tot}^{\epsilon}(t) \geq K) \leq \frac{\sup_{\epsilon} \sup_{t} \mathbb{E}(Y_{tot}^{\epsilon}(t) + Z_{tot}^{\epsilon}(t))}{K} $$

Ainsi avec $B_K = \{(y, z) \in (\mathbb{R}^+)^d \times \mathbb{N}^l, y_{tot} + z_{tot} \leq K\}$, pour tout $\nu > 0$, il existe $K > 0$ tel que

$$ \forall \epsilon, t > 0, \inf_{\epsilon} \mathbb{E}(\Gamma_{\epsilon}([0, t], B_K)) \geq t(1 - \nu) $$

et avec le Lemme 1.3 de [40], nous avons la compacité relative de $(X^\epsilon, \Gamma_\epsilon)$ dans $\mathcal{D}_{\mathbb{R}^+}[0, \infty[ \times \mathcal{L}_m(\mathbb{R}^d \times \mathbb{N}^l)$.

## 4.2 Identification de la limite

Soit $(x, \Gamma)$ un point limite de $\Psi^\epsilon$ dans $\mathcal{D}_{\mathbb{R}^+}[0, \infty[ \times \mathcal{L}_m(\mathbb{R}^d \times \mathbb{N}^l)$.

Dans un premier temps, nous identifions la mesure limite $\Gamma$. Avec $f$ dépendant séparément de $y$ et $z$, c'est-à-dire $f(y, z) = \phi(y)\bar{f}(z)$ :

$$ \begin{aligned} A^{\epsilon} f(x, y, z) = & \frac{\lambda_0(y, z)}{\epsilon} x \bar{f}(z) [\phi(y + \epsilon^p e_1) - \phi(y)] \\ & + \sum_{i=1}^{d} \frac{\mu_i(y, z)}{\epsilon} y_i \bar{f}(z) [\phi(y - \epsilon^p e_i) - \phi(y)] + \sum_{i=1}^{d-1} \frac{\lambda_i(y, z)}{\epsilon} y_i \bar{f}(z) [\phi(y + \epsilon^p (e_{i+1} - e_i)) - \phi(y)] \\ & + \frac{\lambda_d(z)}{\epsilon} y_d [\phi(y - \epsilon^p e_d) \bar{f}(z + e_1) - \phi(y) \bar{f}(z)] \\ & + \sum_{i=1}^{l} \frac{\mu_{i+d}(z)}{\epsilon} z_i \phi(y) [\bar{f}(z - e_i) - \bar{f}(z)] + \sum_{i=1}^{l-1} \frac{\lambda_{i+d}(z)}{\epsilon} z_i \phi(y) [\bar{f}(z + e_{i+1} - e_i) - \bar{f}(z)] \end{aligned} $$

Comme

$$ \phi(y - \epsilon^p e_d) \bar{f}(z + e_1) - \phi(y) \bar{f}(z) = [\phi(y - \epsilon^p e_d) - \phi(y)] \bar{f}(z + e_1) + \phi(y) [\bar{f}(z + e_1) - \bar{f}(z)] $$
avec le générateur $L_{y_d}$ de la partie rapide défini par (2.1), on a :

$$
\begin{align*}
A^\epsilon f(x, y, z) ={}& \frac{\lambda_0(y, z)}{\epsilon} x \bar{f}(z) [\phi(y + \epsilon^p e_1) - \phi(y)] + \sum_{i=1}^{d} \frac{\mu_i(y, z)}{\epsilon} y_i \bar{f}(z) [\phi(y - \epsilon^p e_i) - \phi(y)] \\
& + \sum_{i=1}^{d-1} \frac{\lambda_i(y, z)}{\epsilon} y_i \bar{f}(z) [\phi(y + \epsilon^p (e_{i+1} - e_i)) - \phi(y)] \\
& + \frac{\lambda_d(z)}{\epsilon} y_d [\phi(y - \epsilon^p e_d) - \phi(y)] \bar{f}(z + e_1) \\
& + \frac{1}{\epsilon} \phi(y) L_{y_d} h(z)
\end{align*}
$$

De plus, le processus

$$
\epsilon[f(Y^{\epsilon}(t), Z^{\epsilon}(t)) - f(Y^{\epsilon}(0), Z^{\epsilon}(0))] - \int_{0}^{t} \int_{\mathbb{R}^{d} \times \mathbb{N}^{l}} \epsilon A^{\epsilon} f(X^{\epsilon}(s), y, z) \Gamma_{\epsilon}(ds, dy, dz) \quad (2.7)
$$

est une $\mathcal{F}_t^\epsilon$-martingale.

Avec le Lemme 1.5 de [40], comme $\epsilon A^\epsilon f(X^\epsilon(s), y, z) \xrightarrow{\epsilon \to 0} \phi(y) L_{y_d} \bar{f}(z)$ :

$$
\begin{align*}
\forall \phi C^1, \lim_{\epsilon \to 0} \int_0^t \int_{\mathbb{R}^d \times \mathbb{N}^l} \epsilon A^\epsilon f(X^\epsilon(s), y, z) \Gamma_\epsilon(ds, dy, dz) &= \int_0^t \int_{\mathbb{R}^d \times \mathbb{N}^l} \phi(y) L_{y_d} \bar{f}(z) \Gamma(ds, dy, dz) \\
&= \int_0^t \int_{\mathbb{R}^d} \phi(y) \int_{\mathbb{N}^l} L_{y_d} \bar{f}(z) \Gamma(ds, dy, dz)
\end{align*}
$$

D'après (2.7),

$$
\int_{0}^{t} \int_{\mathbb{R}^{d}} \phi(y) \int_{\mathbb{N}^{l}} L_{y_{d}} \bar{f}(z) \Gamma(ds, dy, dz)
$$

est une martingale continue à variations bornées, par conséquent elle est égale à 0 pour tout $\bar{f}$ et $\phi$.
Avec le Lemme 1.4 de [40], il existe un processus $\gamma_s$ tel que, pour tout $\bar{f}$ et pour tout $\phi$ :

$$
\int_0^t \int_{\mathbb{R}^d} \phi(y) \int_{\mathbb{N}^l} L_{y_d} \bar{f}(z) \Gamma(ds, dy, dz) = \int_0^t \int_{\mathbb{R}^d} \phi(y) \int_{\mathbb{N}^l} L_{y_d} \bar{f}(z) \gamma_s(dy, dz) \, ds
$$

D'après le théorème de désintégration de mesure (Chapitre III de [14]), avec la projection na-
turelle $\pi : \mathbb{R}^d \times \mathbb{R}^l \rightarrow \mathbb{R}^d$, il existe une famille de mesures de probabilité $(\tilde{\gamma}_{s,y})_{y \in \mathbb{R}^d}$ définie
presque partout telle que :

$$
\int_{\mathbb{R}^d} \phi(y) \int_{\mathbb{N}^l} L_{y_d} \bar{f}(z) \gamma_s(dy, dz) = \int_{\mathbb{R}^d} \phi(y) \left[ \int_{\mathbb{N}^l} L_{y_d} \bar{f}(z) \tilde{\gamma}_{s,y}(dz) \right] \tilde{\nu}_s(dy)
$$

Donc nous avons, pour presque tout $s \in [0, T]$,

$$
\forall \phi C^1, \int_{\mathbb{R}^d} \phi(y) \left[ \int_{\mathbb{N}^l} L_{y_d} \bar{f}(z) \tilde{\gamma}_{s,y}(dz) \right] \tilde{\nu}_s(dy) = 0
$$
ce qui implique d'après le théorème de Riesz-Markov-Kakutani (théorème 6.19 de [59]) que
pour presque tout $s \in [0, T]$,

$$
\int_{\mathbb{N}^l} L_{y_d} \bar{f}(z) \tilde{\gamma}_{s,y}(dz) = 0, \tilde{\nu}_s \text{ p.p}
$$

Donc d'après la Proposition 3.2, nous pouvons conclure que pour presque tout $s \in [0, T]$,
$\tilde{\gamma}_{s,y}(dz) = \pi_{y_d}(dz)$, $\tilde{\nu}_s$ p.p.

L'étape suivante consiste à identifier la mesure $\tilde{\nu}_s$. Pour une fonction $h$ dépendant seulement de la seconde variable $y$, c'est-à-dire $\forall (x, z), h(x, \cdot, z) = h(\cdot)$, pour presque tout $s \in [0, T]$ :

$$
\begin{align*}
A^{\epsilon} h(x, y, z) &= \frac{\lambda_0(y, z)}{\epsilon} x[h(y + \epsilon^p e_1) - h(y)] \\
&\quad + \sum_{i=1}^{d} \frac{\mu_i(y, z)}{\epsilon} y_i[h(y - \epsilon^p e_i) - h(y)] \\
&\quad + \sum_{i=1}^{d-1} \frac{\lambda_i(y, z)}{\epsilon} y_i[h(y + \epsilon^p (e_{i+1} - e_i)) - h(y)] \\
&\quad + \frac{\lambda_d(z)}{\epsilon} y_d[h(y - \epsilon^p e_d) - h(y)] \\
&= \frac{1}{\epsilon^{1-p}} \mathcal{H}_x h(x, y, z) + O(1) \tag{2.8}
\end{align*}
$$

avec l'opérateur suivant

$$
\mathcal{H}_x h(x, y, z) = \lambda_0(y, z) x \frac{\partial h}{\partial e_1}(y) - \sum_{i=1}^d \left( \mu_i(y, z) + \lambda_i(y, z) \right) y_i \frac{\partial h}{\partial e_i}(y) + \sum_{i=1}^{d-1} \lambda_i(y, z) y_i \frac{\partial h}{\partial e_{i+1}}(y) \tag{2.9}
$$

Par définition du générateur infinitésimal, le processus :

$$
\epsilon^{1-p} [f(Y^{\epsilon}(t)) - f(Y^{\epsilon}(0))] - \int_{0}^{t} \int_{\mathbb{R}^{d} \times \mathbb{N}^{l}} \mathcal{H}_{X^{\epsilon}(s)} h(X^{\epsilon}(s), y, z) \Gamma_{\epsilon}(ds, dy, dz) + R^{\epsilon} h(t) \quad (2.10)
$$

est une $\mathcal{F}_t^\epsilon$-martingale avec

$$
R^{\epsilon}h(t) = \int_{0}^{t} \left[ \mathcal{H}_{X^{\epsilon}(s)}h(X^{\epsilon}(s), Y^{\epsilon}(s), Z^{\epsilon}(s)) - \epsilon^{1-p}A^{\epsilon}h(X^{\epsilon}(s), Y^{\epsilon}(s), Z^{\epsilon}(s)) \right] ds
$$

En utilisant les mêmes arguments que précédemment

$$
\int_{\mathbb{R}^d \times \mathbb{N}^l} \mathcal{H}_{x(s)} h(x(s), y, z) \tilde{\gamma}_{s,y}(dz) \tilde{\nu}_s(dy) = 0 \text{ presque sûrement}
$$

Alors,

$$
\begin{align*}
\int_{\mathbb{R}^d \times \mathbb{N}^l} \mathcal{H}_{x(s)} h(x(s), y, z) \tilde{\gamma}_{s,y}(dz) \tilde{\nu}_s(dy) &= \int_{\text{Supp}(\tilde{\nu}_s) \times \mathbb{N}^l} \mathcal{H}_{x(s)} h(x(s), y, z) \pi_{y_d}(dz) \tilde{\nu}(dy) \\
&= \int_{\text{Supp}(\tilde{\nu}_s)} \left( \int_{\mathbb{N}^l} \mathcal{H}_{x(s)} h(x(s), y, z) \pi_{y_d}(dz) \right) \tilde{\nu}_s(dy) \\
&= \int_{\mathbb{R}^d} \tilde{\mathcal{H}}_{x(s)} h(x(s), y) \tilde{\nu}_s(dy)
\end{align*}
$$
car $\tilde{\mathcal{H}}_{x(s)}h(\cdot, y)$ est défini sur $\mathbb{R}^d$ (voir **Définition 3.2**).

Avec l'hypothèse 3.2, pour presque tout $s \in [0, T]$,

$\tilde{\nu}_s(dy) = \delta_{v_{x(s)}^*}(dy)$ avec $v^*$ l'unique point fixe défini par la **Proposition 3.3**.

Nous avons donc identifié la mesure limite $\Gamma$ :

$$ \Gamma(ds, dy, dz) = \pi_{(y_{x(s)}^*)(d)}(dz) \delta_{y_{x(s)}^*}(dy) \, ds. $$

Il reste à identifier la fonction $x$. Tout d'abord, comme la première composante $X^\epsilon$ fait des sauts tendant vers 0 (quand $\epsilon$ tend vers 0), c'est-à-dire pour tout $T > 0$,

$$ \sup_{t \in [0,T]} |X^{\epsilon}(t) - X^{\epsilon}(t^{-})| \le \epsilon \to 0, $$

$x$ est une fonction continue.

Soit $A : C_b^2(\mathbb{R}^+) \to C_b^1(\mathbb{R}^+ \times \mathbb{N}^d)$ le générateur infinitésimal défini pour tout $(x, y, z) \in \mathbb{F}$ par :

$$ Af(x, y, z) = (r_0(y, z) - \lambda_0(y, z) - \mu_0(y, z))xf'(x). $$

$$ \text{Avec } R^{\epsilon}(t) = \int_{0}^{t} [AId - A^{\epsilon}Id] (X^{\epsilon}(s), Y^{\epsilon}(s), Z^{\epsilon}(s)) \, ds, $$

$$ M_{Id}^{\epsilon}(t) = f(X^{\epsilon}(t)) - f(X^{\epsilon}(0)) - \int_{0}^{t} AId(X^{\epsilon}(s), Y^{\epsilon}(s), Z^{\epsilon}(s)) \, ds + R^{\epsilon}(t) \quad (2.11) $$

est une $\mathcal{F}_t^\epsilon$-martingale à variation quadratique prévisible

$$ \langle M_{Id}^{\epsilon} \rangle_t = \epsilon \int_0^t [r_0 + \lambda_0 + \mu_0] (Y^{\epsilon}(s), Z^{\epsilon}(s)) \, ds. $$

Avec la **Proposition 3.1**, $\mathbb{E}(\langle M_{Id}^{\epsilon} \rangle_t) \xrightarrow{\epsilon \to 0} 0$ et $M_t \equiv 0$ pour tout $t > 0$.

Alors comme pour le cas à 2 échelles de temps, en appliquant le Lemme 1.5d de [40], il existe une filtration $(\mathcal{F}_t)_t$ telle que pour $f(x, y, z) = x$, le processus

$$ M_t = x(t) - x^{in} - \int_0^t \int_{\mathbb{R}^d \times \mathbb{N}^l} AId(x(s), y, z)\Gamma(ds, dy, dz) \quad (2.12) $$

est une $(\mathcal{F}_t)_t$-martingale. Or pour $s > 0$,

$$ \int_{\mathbb{R}^d \times \mathbb{N}^l} AId(x(s), y, z)\Gamma(ds, dy, dz) = x(s) \int_{\mathbb{R}^d} (R_0(y) - \Lambda_0(y) - M_0(x(s), y)) \delta_{y_{x(s)}^*}(dy) \\ = x(s)D_0(x(s)) $$

Ainsi, $x$ est solution de

$$ x(t) = x^{in} + \int_{0}^{t} x(s)D_{0}(x(s)) \, ds $$

Comme $S_f$ est continue pour toute fonction $f$ bornée par $V_p + 1$ pour $p \in \mathcal{P}$ (**Proposition 3.2**), la fonction $D_0$ est continue et par composition, nous obtenons que $x$ appartient à $\mathcal{C}^1(\mathbb{R}^+)$ et est solution de :

$$ \frac{dx}{dt}(t) = D_0(x(t))x(t), \quad x(0) = x^{in}. $$
## 4.3 Unicité de la limite

L'unicité de la mesure stationnaire $\Gamma$ découle directement de l'unicité de la mesure de probabilité invariante $\pi_{y_d}$ de la partie rapide (Proposition 3.2) et de l'unicité d'une mesure stationnaire pour l'opérateur $\tilde{\mathcal{H}}_x$ (Hypothèse 3.2).

Prouvons l'unicité de $x$, solution de

$$ \frac{dw}{dt} = D_0(w)w, \quad w(0) = x^{in}. \qquad (2.13) $$

Montrons que la fonction $F$ définie sur $\mathbb{R}^+$ par

$$ F(w) = w\langle r_0(y_w^*, \cdot) - \lambda_0(y_w^*, \cdot) - \mu_0(y_w^*, \cdot), \pi_{(y_w^*)d} \rangle $$

est localement lipschitzienne.

Soient $w_1 \in \mathbb{R}^+$ et $U \subset \mathbb{R}^+$ un voisinage de $w_1$. Pour tout $w_2 \in U$, avec $g(w, \cdot) = r_0(y_w^*, \cdot) - \lambda_0(y_w^*, \cdot) - \mu_0(y_w^*, \cdot)$,

$$ F(w_1) - F(w_2) = w_1\langle g(w_1, \cdot), \pi_{(y_{w_1}^*)d} - \pi_{(y_{w_2}^*)d} \rangle + (w_1 - w_2)\langle g(w_1, \cdot), \pi_{(y_{w_2}^*)d} \rangle + w_2\langle g(w_1, \cdot) - g(w_2, \cdot), \pi_{(y_{w_2}^*)d} \rangle \quad (2.14) $$

Avec les Hypothèses 2.3 et 2.6, il existe $p \in \mathcal{P}$ et $\alpha_{w_1} > 0$ tels que $g(w_1, \cdot) \le \alpha_{w_1}(1 + V_p)$. Avec le caractère lipschitzien de la fonction $S_g$, il reste seulement à traiter le dernier terme de (2.14) pour pouvoir conclure que la fonction $F$ est localement lipschitzienne.

Or avec les Hypothèses 2.6 et 3.3, la fonction $g$ est localement lipschitzienne selon sa première variable, donc il existe $C > 0$ (dépendant de $w_1$) tel que

$$ \forall w_2 \in U, \quad \left|F(w_1) - F(w_2)\right| \le C|w_1 - w_2| $$

D'après le théorème de Cauchy-Lipschitz, la solution maximale de l'équation (2.13) est unique. Enfin, comme la fonction $D_0$ est bornée supérieurement par $R_0$ (Hypothèse 2.3), il y a une unique solution globale à (2.13).

# 5 Discussion du modèle limite

Les outils mathématiques utilisés afin de prouver l'existence et l'unicité d'une limite au modèle à trois échelles de temps sont les mêmes que ceux du chapitre 1. Néanmoins, nous avons dû faire les hypothèses supplémentaires suivantes :

* une hypothèse sur la dépendance des taux de transition permettant de séparer plus facilement chaque échelle de temps (Hypothèse 2.2),

* une hypothèse sur la monotonie des taux intervenant au sein de la partie intermédiaire du processus (Hypothèse 3.1) afin de garantir l'existence d'un unique point fixe au système (2.2),

* une hypothèse sur la mesure stationnaire de l'opérateur $\tilde{\mathcal{H}}_x$ (Hypothèse 3.2).

* une hypothèse sur la dépendance du point fixe $y^*$ à l'effectif de la réserve (Hypothèse 3.3).
La première hypothèse est nécessaire à la fois pour que chaque échelle de temps interagisse simplement avec les autres échelles et pour être capable d'étudier la distribution limite $\pi$ et le point fixe $y^*$. D'un point de vue biologique, il n'est pas choquant qu'un follicule régule seulement les follicules dans le même stade de maturation (effet de compétition) et les follicules moins matures (ralentissement de la croissance pour préserver la capacité de reproduction le plus longtemps possible).

Les Hypothèses 3.1 et 3.2 sont étroitement liées à l'étude du système d'EDOs caractérisant l'échelle de temps intermédiaire du modèle. En dimension quelconque, les points fixes de tels systèmes sont particulièrement ardus à déterminer ce qui nous oblige à faire des hypothèses aussi contraignantes. En très petite dimension ($d$ égal 1 ou 2) ou avec des formes de taux de transition connues, l'hypothèse sur la monotonie des taux de l'échelle intermédiaire pourrait probablement être relaxée.

Concernant la mesure stationnaire de l'opérateur $\tilde{\mathcal{H}}_x$, la situation est un peu différente. Deux propriétés doivent être vérifiées : la non-existence d'un cycle limite et l'unicité d'une mesure stationnaire pour l'opérateur $\tilde{\mathcal{H}}_x$.

L'hypothèse suivante garantit que le point fixe $y^*$ est atteint quelque soit la condition initiale, excluant ainsi l'existence d'un cycle limite :

**Hypothèse 5.1**

$\forall x > 0$, pour toute condition initiale $v_0 \in (\mathbb{R}^+)^d$, $\|v(t) - v^*_x\| \underset{t \to \infty}{\longrightarrow} 0$

avec $v(t)$ la solution du système (2.2).

La principale difficulté pour pouvoir utiliser cette hypothèse est la traduction d'une propriété sur le système d'EDO en une propriété sur les mesures stationnaires de l'opérateur associé à ce système. Ce travail est encore en cours afin de remplacer l'hypothèse 3.2 par l'hypothèse 5.1.

Enfin la dernière hypothèse sur la dépendance (en $x$) du point fixe $y^*$ est une hypothèse technique. Cette hypothèse est cohérente avec le fait que l'opérateur associé à ce point fixe est linéaire selon la variable $x$. Malgré la simplicité de cette opérateur et l'hypothèse faite sur la monotonie des taux, nous n'avons pas réussi à prouver la dépendance linéaire en $x$ du point fixe $y^*$, ce qui nous a conduit à supposer ce résultat.

# 6 Étude de l'erreur d'approximation

Dans la suite, nous nous intéressons à l'erreur entre le modèle et sa limite dans les cas où l'on considère 2 (Chapitre 1) ou 3 (Chapitre 2) échelles de temps. La méthode utilisée est inspirée de [36] et ce travail a été effectué en collaboration avec Léa Popovic<fnref n="1" />.

Afin de présenter la méthode utilisée, nous présentons tout d'abord le cas avec deux échelles de temps.

## 6.1 Erreur d'approximation pour deux échelles de temps

Nous utilisons les mêmes notations que celles du chapitre 1. Sous certaines hypothèses, le processus $(X^\epsilon, Y^\epsilon)$ converge vers le processus $(x, \Gamma)$ où $x$ est solution d'une équation différentielle

<footnote n="1">1. Department of Mathematics & Statistics, Concordia University, Montréal, Canada</footnote>
ordinaire couplée avec la mesure $\Gamma$.

La convergence de la seconde composante $Y^\epsilon$ a lieu à travers la convergence de sa mesure d'occupation vers $\Gamma$. Nous étudierons donc seulement l'erreur entre la réserve $X^\epsilon$ et sa limite $x$.

### 6.1.1 Décomposition de l'erreur

Soit $U^\epsilon(t) = \frac{X^\epsilon(t)-x(t)}{\sqrt{\epsilon}}$ l'erreur entre la réserve et sa limite.

Soient $F(x, y) = f(y)x$ le taux de décroissance de la réserve avec $f(y) = r_0(y) - \lambda_0(y) - \mu_0(y)$.

D'après la définition du générateur infinitésimal, le processus suivant :

$$M^{\epsilon,1}(t) = X^{\epsilon}(t) - X^{\epsilon}(0) - \int_{0}^{t} F^{\epsilon}(X^{\epsilon}(s), Y^{\epsilon}(s)) \, ds, \quad (2.15)$$

est une $\mathcal{F}_t^\epsilon$-martingale.

Pour une fonction $H^\epsilon$ bien choisie, l'erreur $U^\epsilon$ peut se décomposer de la façon suivante :

$$\begin{align}
U^{\epsilon}(t) &= U^{\epsilon}(0) + \frac{1}{\sqrt{\epsilon}} \int_{0}^{t} F(X^{\epsilon}(s), Y^{\epsilon}(s)) \, ds - \frac{1}{\sqrt{\epsilon}} \int_{0}^{t} \Lambda_{0}(x(s)) x(s) \, ds + \frac{1}{\sqrt{\epsilon}} M^{\epsilon,1}(t) \nonumber \\
&= U^{\epsilon}(0) + \frac{1}{\sqrt{\epsilon}} \int_{0}^{t} \left[ \Lambda_{0}(X^{\epsilon}(s)) X^{\epsilon}(s) - \Lambda_{0}(x(s)) x(s) \right] ds \tag{2.16}
\end{align}$$

$$+ \frac{1}{\sqrt{\epsilon}} \int_{0}^{t} \left[ F(X^{\epsilon}(s), Y^{\epsilon}(s)) - \Lambda_{0}(X^{\epsilon}(s)) X^{\epsilon}(s) - A^{\epsilon} H^{\epsilon}(X^{\epsilon}(s), Y^{\epsilon}(s)) \right] ds \quad (2.17)$$

$$+ \frac{1}{\sqrt{\epsilon}} \left[ H^{\epsilon}(X^{\epsilon}(t), Y^{\epsilon}(t)) - H^{\epsilon}(X^{\epsilon}(0), Y^{\epsilon}(0)) \right] \quad (2.18)$$

$$+ \frac{1}{\sqrt{\epsilon}} (M^{\epsilon,1}(t) - M^{\epsilon,2}(t)) \quad (2.19)$$

avec la $\mathcal{F}_t^\epsilon$-martingale suivante :

$$M^{\epsilon,2}(t) = H^{\epsilon}(X^{\epsilon}(t), Y^{\epsilon}(t)) - H^{\epsilon}(X^{\epsilon}(0), Y^{\epsilon}(0)) - \int_{0}^{t} A^{\epsilon} H^{\epsilon}(X^{\epsilon}(s), Y^{\epsilon}(s)) \, ds. \quad (2.20)$$

La caractérisation de l'erreur $U^\epsilon$ se fait en trois étapes : le choix d'une fonction $H^\epsilon$ afin que les termes (2.17) and (2.18) tendent vers 0 quand $\epsilon$ tend vers 0, puis le calcul de la variation quadratique de la différence des martingales (terme (2.19)) et enfin le passage à la limite dans l'expression ci-dessus.

### 6.1.2 Sélection de la fonction $H^\epsilon$

D'après la définition de $\Lambda_0$, pour tout $(x, y) \in \mathbb{R}^+ \times (\mathbb{R}^+)^d$,

$$F(x, y) - \Lambda_0(x)x = x \left[ (r_0(y) - \lambda_0(y) - \mu_0(y)) - \langle r_0 - \lambda_0 - \mu_0, \pi_x \rangle \right]$$

Avec $f = r_0 - \lambda_0 - \mu_0$, cette expression est semblable à l'équation de Poisson associée à $(f, L_x)$ :

$$L_x g(y) = \langle f, \pi_x \rangle - f(y) \quad (2.21)$$
D'après la preuve du théorème 3.1 du chapitre 1, pour tout $x > 0$, l'équation de Poisson admet une
solution $g_x$ tel qu'il existe $\alpha_x < \infty$ tel que $|g_x| \le \alpha_x(V_p + 1)$.

Soient $\tilde{g}_x = g_x(\cdot + e_1)$ et $\Delta_x = (\tilde{g}_x - g_x)\lambda_0$. Après calculs, nous vérifions que $\frac{g_{x+\varepsilon} - g_x}{\varepsilon}$ est solution de l'équation de Poisson associée à $(\Delta_x, L_{x+\varepsilon})$.

Ainsi, pour tout compact $K$ de $\mathbb{N}^d$,

$$
\forall x > 0, \forall \varepsilon > 0, \forall y \in K, g_{x+\varepsilon}(y) - g_x(y) = \mathcal{O}(\varepsilon) \quad (2.22)
$$

La proposition suivante indique le choix de la fonction $H^\varepsilon$ pour que le terme (2.17) dans la décomposition de l'erreur tende vers 0.

**Proposition 6.1**

Soient pour tout $(x, y) \in \mathbb{R}^{+} \times \mathbb{N}^{d}$, $H^{\epsilon}(x, y) = -\epsilon x g_{x}(y)$ et

$$
I^{\epsilon}(x, y) = F^{\epsilon}(x, y) - \Lambda_{0}(x)x - A^{\epsilon} H^{\epsilon}(x, y)
$$

Alors, pour tout compact $K \subset \mathbb{R}^{+} \times \mathbb{N}^{d}$,

$$
\lim_{\epsilon \to 0} \sup_{(x,y) \in K} \frac{1}{\sqrt{\epsilon}} |I^{\epsilon}(x,y)| = 0
$$

**Démonstration :** Soient $K$ un compact de $\mathbb{R}^+ \times \mathbb{N}^d$ et $(x, y)$ un élément de ce compact. D'après la
définition du générateur $A^\epsilon$ (Definition 2.1 du chapitre 1) :

$$
\begin{align*}
A^{\epsilon} H^{\epsilon}(x, y) &= r_0(y) x \left[ x g_x(y) - (x + \epsilon) g_{x+\epsilon}(y) \right] + \lambda_0(y) x \left[ x g_x(y) - (x - \epsilon) g_{x-\epsilon}(y + e_1) \right] \\
&\quad + \mu_0(y) x \left[ x g_x(y) - (x - \epsilon) g_{x-\epsilon}(y) \right] + \sum_{i=1}^{d} \mu_i(y) y_i x \left[ g_x(y) - g_x(y - e_i) \right] \\
&\quad + \sum_{i=1}^{d-1} \lambda_i(y) y_i x \left[ g_x(y) - g_x(y + e_{i+1} - e_i) \right]. \\
&= r_0(y) x^2 \left[ g_x(y) - g_{x+\epsilon}(y) \right] + \lambda_0(y) x^2 \left[ g_x(y) - g_{x-\epsilon}(y + e_1) \right] \\
&\quad + \mu_0(y) x^2 \left[ g_x(y) - g_{x-\epsilon}(y) \right] - x \left[ L_x g_x(y) - \lambda_0(y) x (g_x(y + e_1) - g_x(y)) \right] \\
&\quad + \epsilon \left[ -r_0(y) x g_{x+\epsilon}(y) + \lambda_0(y) x g_{x-\epsilon}(y + e_1) + \mu_0(y) x g_{x-\epsilon}(y) \right]
\end{align*}
$$

Donc comme $g_x$ est solution de l'équation de Poisson associée à $(f, L_x)$ :

$$
I^{\epsilon}(x, y) = -r_{0}(y)x^{2}[g_{x}(y) - g_{x+\epsilon}(y)] - \lambda_{0}(y)x^{2}[g_{x}(y + e_{1}) - g_{x-\epsilon}(y + e_{1})] \\
-\mu_{0}(y)x^{2}[g_{x}(y) - g_{x-\epsilon}(y)] + \epsilon[r_{0}(y)xg_{x+\epsilon}(y) - \lambda_{0}(y)xg_{x-\epsilon}(y + e_{1}) - \mu_{0}(y)xg_{x-\epsilon}(y)]
$$

Avec l'équation (2.22), les termes oranges tendent vers 0 à vitesse $\epsilon$. De plus, comme $K$ est un
compact et que toutes les fonctions apparaissant dans $I^\epsilon$ sont bornées sur un compact :

$$
\lim_{\epsilon \to 0} \sup_{x,y} \frac{1}{\sqrt{\epsilon}} |I^{\epsilon}(x,y)| = 0
$$

<figure><img src="image_26.png" /></figure>
### 6.1.3 Variation quadratique

Soit $M^\epsilon = \frac{1}{\sqrt{\epsilon}}(M^{\epsilon,1} - M^{\epsilon,2})$.

Afin de caractériser la variance de la limite de la martingale $M^\epsilon$, nous devons identifier la fonction $G$ telle que

$$ \frac{1}{\epsilon} [X^{\epsilon} - H^{\epsilon}(X^{\epsilon}, Y^{\epsilon})]_{t} - \int_{0}^{t} G(X^{\epsilon}(s), Y^{\epsilon}(s)) \, ds \to 0 \quad (2.23) $$

Pour cela, nous calculons la variation quadratique de $X^\epsilon - H^\epsilon(X^\epsilon, Y^\epsilon)$, la multiplions par $\frac{1}{\epsilon}$ et identifions les termes d'ordre 1.

Avec l'expression précédente de $H^\epsilon(X^\epsilon, Y^\epsilon)$, nous avons

$$ \begin{aligned} \frac{1}{\epsilon} [X^{\epsilon} - H^{\epsilon}(X^{\epsilon}, Y^{\epsilon})]_{t} &= \frac{1}{\epsilon} [X^{\epsilon} + \epsilon X^{\epsilon} g_{X^{\epsilon}}(Y^{\epsilon})]_{t} \\ &= \int_{0}^{t} [\epsilon + \epsilon (X^{\epsilon} + \epsilon) g_{X^{\epsilon}+\epsilon}(Y^{\epsilon}) - \epsilon X^{\epsilon} g_{X^{\epsilon}}(Y^{\epsilon})]^{2} \frac{r_{0}(Y^{\epsilon})}{\epsilon^{2}} X^{\epsilon} \\ &\quad + \int_{0}^{t} [-\epsilon + \epsilon (X^{\epsilon} - \epsilon) g_{X^{\epsilon}-\epsilon}(Y^{\epsilon}) - \epsilon X^{\epsilon} g_{X^{\epsilon}}(Y^{\epsilon})]^{2} \frac{\mu_{0}(Y^{\epsilon})}{\epsilon^{2}} X^{\epsilon} \\ &\quad + \int_{0}^{t} [-\epsilon + \epsilon (X^{\epsilon} - \epsilon) g_{X^{\epsilon}-\epsilon}(Y^{\epsilon} + e_{1}) - \epsilon X^{\epsilon} g_{X^{\epsilon}}(Y^{\epsilon})]^{2} \frac{\lambda_{0}(Y^{\epsilon})}{\epsilon^{2}} X^{\epsilon} \\ &\quad + \sum_{i=1}^{d} \int_{0}^{t} [\epsilon X^{\epsilon} (g_{X^{\epsilon}}(Y^{\epsilon} - e_{i}) - g_{X^{\epsilon}}(Y^{\epsilon}))]^{2} \frac{\mu_{i}(Y^{\epsilon})}{\epsilon^{2}} Y_{i}^{\epsilon} \\ &\quad + \sum_{i=1}^{d-1} \int_{0}^{t} [\epsilon X^{\epsilon} (g_{X^{\epsilon}}(Y^{\epsilon} + e_{i+1} - e_{i}) - g_{X^{\epsilon}}(Y^{\epsilon}))]^{2} \frac{\lambda_{i}(Y^{\epsilon})}{\epsilon^{2}} Y_{i}^{\epsilon} \end{aligned} $$

Comme tous les crochets sont mis au carré, les termes d'ordre 1 sont

$$ \begin{aligned} \frac{1}{\epsilon} [X^{\epsilon} - H^{\epsilon}(X^{\epsilon}, Y^{\epsilon})]_{t} &= \int_{0}^{t} r_{0}(Y^{\epsilon}) X^{\epsilon} + \int_{0}^{t} \mu_{0}(Y^{\epsilon}) X^{\epsilon} \\ &\quad + \int_{0}^{t} [-1 + X^{\epsilon} (g_{X^{\epsilon}}(Y^{\epsilon} + e_{1}) - g_{X^{\epsilon}}(Y^{\epsilon}))]^{2} \lambda_{0}(Y^{\epsilon}) X^{\epsilon} \\ &\quad + \sum_{i=1}^{d} \int_{0}^{t} [X^{\epsilon} (g_{X^{\epsilon}}(Y^{\epsilon} - e_{i}) - g_{X^{\epsilon}}(Y^{\epsilon}))]^{2} \mu_{i}(Y^{\epsilon}) Y_{i}^{\epsilon} \\ &\quad + \sum_{i=1}^{d-1} \int_{0}^{t} [X^{\epsilon} (g_{X^{\epsilon}}(Y^{\epsilon} + e_{i+1} - e_{i}) - g_{X^{\epsilon}}(Y^{\epsilon}))]^{2} \lambda_{i}(Y^{\epsilon}) Y_{i}^{\epsilon} + o(1) \end{aligned} $$

Ainsi pour tout $x, y$, la fonction $G$ solution de (2.23) est :

$$ G(x, y) = r_0(y)x + \mu_0(y)x + [x(g_x(y + e_1) - g_x(y)) - 1]^2 \lambda_0(y)x \quad (2.24) $$

$$ + \sum_{i=1}^{d} [x(g_x(y - e_i) - g_x(y))]^2 \mu_i(y)y_i + \sum_{i=1}^{d-1} [x(g_x(y + e_{i+1} - e_i) - g_x(y))]^2 \lambda_i(y)y_i \quad (2.25) $$
### 6.1.4 Passage à la limite

Avec les éléments des sous-sections précédentes, nous pouvons énoncer le théorème caractérisant l'erreur :

**Théorème 6.1**

Pour tout $t > 0$, $M^\varepsilon(t) = \frac{1}{\sqrt{\varepsilon}}(M^{\varepsilon,1}(t) - M^{\varepsilon,2}(t)) \Rightarrow M$ avec $M$ un processus gaussien de moyenne nulle avec des incréments indépendants :

$$ \mathbb{E}(M(t)M^T(t)) = \int_0^t \langle G(x(s), \cdot), \pi_{x(s)} \rangle \mathrm{d}s $$

Pour toute condition initiale $U^\varepsilon(0)$ telle que, pour $U^{in} \in \mathbb{R}^+$,

$$ U^\varepsilon(0) \underset{a.s}{\rightarrow} U^{in} $$

le processus $U^\varepsilon$ converges vers le processus $U$ vérifiant

$$ U(t) = U^{in} + \int_{0}^{t} \left[ \Lambda'_{0}(x(s))x(s) + \Lambda_{0}(x(s)) \right] U(s) \mathrm{d}s + \int_{0}^{t} \sqrt{\langle G(x(s), \cdot), \pi_{x(s)} \rangle} \mathrm{d}W_{s} $$

avec $W$ le mouvement brownien standard.

**Démonstration :** La martingale $M^\varepsilon$ s'écrit :

$$ M^{\varepsilon}(t) = \frac{1}{\sqrt{\varepsilon}} [X^{\varepsilon}(t) - H^{\varepsilon}(X^{\varepsilon}(t), Y^{\varepsilon}(t))] - \frac{1}{\sqrt{\varepsilon}} [X^{\varepsilon}(0) - H^{\varepsilon}(X^{\varepsilon}(0), Y^{\varepsilon}(0))] \\ - \frac{1}{\sqrt{\varepsilon}} \int_{0}^{t} [A^{\varepsilon} H^{\varepsilon}(X^{\varepsilon}(s), Y^{\varepsilon}(s)) - F(X^{\varepsilon}(s), Y^{\varepsilon}(s)) \mathrm{d}s] $$

Or pour $s > 0$,

$$ M^{\varepsilon}(s) - M^{\varepsilon}(s-) = \frac{1}{\sqrt{\varepsilon}} [X^{\varepsilon}(s) - X^{\varepsilon}(s-)] - \frac{1}{\sqrt{\varepsilon}} [H^{\varepsilon}(X^{\varepsilon}(s), Y^{\varepsilon}(s)) - H^{\varepsilon}(X^{\varepsilon}(s-), Y^{\varepsilon}(s-))] $$

Comme $X^\varepsilon(s) \in \varepsilon\mathbb{N}$ et comme $H^\varepsilon$ est d'ordre $\varepsilon$,

$$ \lim_{\varepsilon \to 0} \mathbb{E} \left[ \sup_{s \le t} | \tilde{M}^{\varepsilon}(s) - \tilde{M}^{\varepsilon}(s-) | \right] = 0 $$

Avec l'expression de la variation quadratique de $M^\varepsilon$ obtenue dans la partie précédente, d'après le théorème central limite pour les martingales (Théorème 7.1.4 de [21]), la martingale $M^\varepsilon$ converges vers le processus Gaussien $M$ avec

$$ \begin{align*} \mathbb{E}(M(t)M^T(t)) &= \lim_{\varepsilon \to 0} \int_0^t G(X^\varepsilon(s), Y^\varepsilon(s)) \mathrm{d}s \\ &= \int_0^t G(x(s), y)\Gamma(dy, \mathrm{d}s) \\ &= \int_0^t \langle G(x(s), \cdot), \pi_{x(s)} \rangle \mathrm{d}s \end{align*} $$

En reprenant la décomposition de $U^\varepsilon$, d'après la définition de $H^\varepsilon$ et avec la Proposition 6.1, les termes (2.17) and (2.18) tendent vers 0 quand $\varepsilon$ tend vers 0.
Comme $X^\epsilon = x + \sqrt{\epsilon}U^\epsilon$, avec un développement de Taylor,

$$ \Lambda_0(X^\epsilon(s))X^\epsilon(s) - \Lambda_0(x(s))x(s) = \sqrt{\epsilon}[\Lambda'_0(x(s))x(s) + \Lambda_0(x(s))]U^\epsilon(s) + \mathcal{O}(\epsilon) $$

Ainsi le terme intégral dans (2.16) converges vers

$$ \int_0^t [\Lambda'_0(x(s))x(s) + \Lambda_0(x(s))] U(s) ds $$

Nous obtenons ainsi que la limite $U$ du processus $U^\epsilon$ est définie par

$$ U(t) = U^{in} + \int_{0}^{t} [\Lambda'_{0}(x(s))x(s) + \Lambda_{0}(x(s))] U(s) ds + \int_{0}^{t} \sqrt{\langle G(x(s), \cdot), \pi_{x(s)} \rangle} dW_{s} $$

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

## 6.2 Erreur d'approximation pour trois échelles de temps

Cette dernière section présente des travaux en cours non finalisé sur l'erreur entre le modèle à trois échelles de temps et sa limite, utilisant la même méthode que précédemment. Rappelons le résultat obtenu dans la section 4 :

### Théorème 6.2

Pour toute condition initiale $\Psi^\epsilon(0) \in \mathbb{F}_0$ (avec $x^{in}$ la limite de $X^\epsilon(0)$), le processus $(X^\epsilon, \Gamma_\epsilon)$ est relativement compact dans $\mathcal{D}_{\mathbb{R}^+}[0, \infty[\times\mathcal{L}_m(\mathbb{R}^d \times \mathbb{N}^l))$, et a un unique processus limite $(x, \Gamma) \in \mathcal{C}^1(\mathbb{R}^+) \times \mathcal{L}_m(\mathbb{R}^d \times \mathbb{N}^l)$ tel que :

$$ \left\{ \begin{aligned} x'(t) &= D_0(x(t))x(t), \quad t > 0, \quad \text{with } x(0) = x^{in} \\ D_0(x) &= \sum_{z \in \mathbb{N}^l} (r_0(y_x^*, z) - \lambda_0(y_x^*, z) - \mu_0(y_x^*, z)) \pi_{(y_x^*)_d}(z) \\ \Gamma(dy, dz, ds) &= \delta_{y_{x(s)}^*}(dy) \pi_{(y_{x(s)}^*)_d}(dz) ds \end{aligned} \right. $$

où,

* pour tout $x > 0$, $y_x^*$ est l'unique point fixe du système d'équations (2.2),

* pour $y \in \mathbb{R}^d$, $\pi_y$ est l'unique mesure de probabilité invariante du processus défini par le générateur $L_{y_d}$ (Proposition 3.2).

L'opérateur associé au système d'équations (2.2) est le suivant, pour $x > 0$ fixé, $h \in \mathcal{C}^1(\mathbb{R}^d)$ :

$$ \tilde{\mathcal{H}}_x h(y) = \Lambda_0(y) x \frac{\partial h}{\partial e_1}(y) - \sum_{i=1}^{d} (M_i(y) + \Lambda_i(y)) y_i \frac{\partial h}{\partial e_i}(y) + \sum_{i=1}^{d-1} \Lambda_i(y) y_i \frac{\partial h}{\partial e_{i+1}}(y) $$

Soit $U^\epsilon(t) = \frac{X^\epsilon(t)-x(t)}{\sqrt{\epsilon}}$ l'erreur entre la réserve $X^\epsilon$ et sa limite $x$. Soient les trois fonctions suivantes :

$$ \left\{ \begin{aligned} f(y, z) &= r_0(y, z) - \lambda_0(y, z) - \mu_0(y, z) \\ F(x, y, z) &= x f(y, z) \\ \bar{F}(x, y) &= x \langle f(y, \cdot), \pi_y \rangle \end{aligned} \right. $$

La caractérisation de l'erreur d'approximation (c'est-à-dire la limite de $U^\epsilon$ quand $\epsilon$ tend vers 0) est obtenue de la façon suivante :
* le processus $U^\epsilon$ est décomposé en 5 termes en introduisant deux fonction $H_1^\epsilon$ et $H_2^\epsilon$,

* les fonctions $H_1^\epsilon$ et $H_2^\epsilon$ sont choisies à l'aide de deux équations de Poisson et une équation aux dérivées partielles afin que trois termes parmi les 5 tendent vers 0,

* le théorème central limite pour les martingales permet de conclure sur la limite de $U^\epsilon$.

Dans toute la suite, nous supposerons que la **vitesse de l'échelle intermédiaire** $\epsilon^p$ est telle que $p > \frac{1}{2}$. Cette hypothèse est discutée dans la sous-section 6.3.

### 6.2.1 Décomposition de l'erreur

Comme $\bar{F}(x, y_x^*) = xD_0(x)$,

$$U^{\epsilon}(t) = U^{\epsilon}(0) + \frac{1}{\sqrt{\epsilon}} \int_{0}^{t} \left[ F(\Psi^{\epsilon}(s)) - \bar{F}(x(s), y_{x(s)}^{*}) \right] ds + \frac{1}{\sqrt{\epsilon}} M^{\epsilon,1}(t)$$

où

$$M^{\epsilon,1}(t) = X^{\epsilon}(t) - X^{\epsilon}(0) - \int_{0}^{t} F(\Psi^{\epsilon}(s)) ds, \qquad (2.26)$$

est une $\mathcal{F}_t^\epsilon$-martingale.

Pour deux fonctions à déterminer $H_1^\epsilon$ et $H_2^\epsilon$, nous décomposons l'erreur $U^\epsilon$ de la façon suivante :

$$U^{\epsilon}(t) = U^{\epsilon}(0) + \frac{1}{\sqrt{\epsilon}} \int_{0}^{t} \left[ \bar{F}\left(X^{\epsilon}(s), y_{X^{\epsilon}(s)}^{*}\right) - \bar{F}\left(x(s), y_{x(s)}^{*}\right) \right] ds \qquad (2.27)$$

$$+ \frac{1}{\sqrt{\epsilon}} \int_{0}^{t} \left[ F(\Psi^{\epsilon}(s)) - \bar{F}\left(X^{\epsilon}(s), Y^{\epsilon}(s)\right) - A^{\epsilon} H_{1}^{\epsilon}(\Psi^{\epsilon}(s)) \right] ds \qquad (2.28)$$

$$+ \frac{1}{\sqrt{\epsilon}} \int_{0}^{t} \left[ \bar{F}\left(X^{\epsilon}(s), Y^{\epsilon}(s)\right) - \bar{F}\left(X^{\epsilon}(s), y_{X^{\epsilon}(s)}^{*}\right) - A^{\epsilon} H_{2}^{\epsilon}(\Psi^{\epsilon}(s)) \right] ds \qquad (2.29)$$

$$+ \frac{1}{\sqrt{\epsilon}} \left[ H_{\text{tot}}^{\epsilon}(\Psi^{\epsilon}(t)) - H^{\text{tot}}(\Psi^{\epsilon}(0)) \right] \qquad (2.30)$$

$$+ \frac{1}{\sqrt{\epsilon}} \left( M^{\epsilon,1}(t) - M^{\epsilon,2}(t) \right) \qquad (2.31)$$

avec $H_{\text{tot}}^\epsilon = H_1^\epsilon + H_2^\epsilon$ et $M^{\epsilon,2}$ la martingale suivante :

$$M^{\epsilon,2}(t) = H_{\text{tot}}^{\epsilon}(\Psi^{\epsilon}(t)) - H_{\text{tot}}^{\epsilon}(\Psi^{\epsilon}(0)) - \int_{0}^{t} A^{\epsilon} H_{\text{tot}}^{\epsilon}(\Psi^{\epsilon}(s)) ds$$

Dans la décomposition ci-dessus, chaque source d'erreur a été séparée et peut ainsi être traitée séparément. Chacune de ses erreurs provient de la différence au sein d'une échelle de temps entre le modèle initial et le modèle limite.

Le terme (2.28) représente l'effet sur la réserve de la différence entre le taux de décroissance effectif de $X^\epsilon(F(\Psi^\epsilon))$ et le taux de décroissance moyennisé selon $\pi(\bar{F}(X^\epsilon, Y^\epsilon))$.

L'erreur provenant de l'échelle intermédiaire ; c'est-à-dire entre le taux de décroissance moyennisé selon $\pi(\bar{F}(X^\epsilon, Y^\epsilon))$ et le même taux évalué au point fixe $y^*$; est représentée par le terme (2.29).

Enfin le terme (2.27) représente l'erreur due à la réserve elle-même.

Les termes (2.30) and (2.31) sont des termes de compensation.