### 6.2.2 Contrôle de l'erreur issue de l'échelle rapide

Cette erreur provient de la différence entre $Z^\epsilon$ et sa limite $\pi_{Y_d^\epsilon}$. La fonction $H_1^\epsilon$ recherchée pour que le terme (2.28) tende vers 0 dépend ainsi seulement de la dernière composante du modèle $z$. De plus, pour tout $(x, y, z) \in \mathbb{F}$, on a

$$F(x, y, z) - \bar{F}(x, y) = x(f(y, z) - \langle f(y, \cdot), \pi_y \rangle)$$

ce qui nous incite à définir l'équation de Poisson associée à $(L_{y_d}, f)$, pour tout $y \in (\mathbb{R}_*^+)^d$ :

$$L_{y_d} g(z) = \langle f(y, \cdot), \pi_y \rangle - f(y, z) \quad (2.32)$$

et pour tout $p \in \mathcal{P}$, la fonction $V_p$ définie pour tout $z \in \mathbb{N}^l$ par

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

Dans la suite, nous aurons besoin d'étudier la même équation de Poisson avec une fonction $f$ différente. Pour cette raison, la proposition suivante est vraie pour toute fonction bornée par $V_p+1$ :

**Proposition 6.2**

Soit une fonction $f$ telle qu'il existe $c > 0$ et $p \in \mathcal{P}$ tels que

$$\forall y \in \mathbb{R}^d, |f(y, \cdot)| \le c(V_p + 1)$$

Alors l'équation de Poisson associée à $(L_{y_d}, f)$ a une solution $g_y$ telle qu'il existe $\alpha_y < \infty$ tel que $|g_y| \le \alpha_y(V_p + 1)$.

**Démonstration:** En utilisant le même raisonnement que dans la démonstration du théorème 3.1 du chapitre 1, on obtient le résultat voulu. ■

La fonction $f$ de la proposition ci-dessus peut dépendre de $y$. Dans ce cas, en s'inspirant du résultat démontré dans le cas des 2 échelles de temps et comme le générateur $L_{y_d}$ est linéaire en $y_d$, nous faisons l'hypothèse suivante (voir discussion en sous-section 6.3) :

**Hypothèse 6.1 (Dépendance en $y$ de la solution de l'équation de Poisson)**

Pour tout compact $K$ de $(\mathbb{R}^+)^d \times \mathbb{N}^l$ et pour tout $\varepsilon > 0$, nous supposons que

$$\forall (y, z) \in K, g_{y+\varepsilon}(z) - g_y(z) = \mathcal{O}(\varepsilon)$$

A l'aide de la solution de l'équation de Poisson, nous choisissons $H_1^\epsilon(x, y, z) = -\varepsilon x g_y(z)$, et en adaptant les calculs de la Proposition 6.1, nous en déduisons la proposition suivante :

**Proposition 6.3**

Soient $g_y$ la solution de l'équation de Poisson associée à $(L_y, f(y, \cdot))$ et pour tout $(x, y, z) \in \mathbb{F}$, $H_1^\epsilon(x, y, z) = -\varepsilon x g_y(z)$.

Alors, pour tout compact $K \subset \mathbb{F}$,

$$\lim_{\varepsilon \to 0} \sup_{(x,y,z) \in K} \frac{1}{\sqrt{\varepsilon}} |F(x, y, z) - \bar{F}(x, y) - A^\varepsilon H_1^\varepsilon(x, y, z)| = 0$$
### 6.2.3 Contrôle de l'erreur issue de l'échelle intermédiaire

Le but de cette partie est de trouver une fonction $H_2^\varepsilon$ telle que le terme (2.29) tende vers 0. Celui-ci provient de la différence entre $Y^\varepsilon$ et $y_{X^\varepsilon}^*$ et est donc lié à l'opérateur $\tilde{\mathcal{H}}_x$ (2.2) associé au système d'EDOs permettant de définir le point fixe $y_{X^\varepsilon}^*$.

**Hypothèse 6.2**

Rappelons que :

$$ \tilde{\mathcal{H}}_x h(y) = \Lambda_0(y) x \frac{\partial h}{\partial e_1}(y) - \sum_{i=1}^{d} \left( M_i(y) + \Lambda_i(y) \right) 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) $$

Soient $x > 0$ et l'équation aux dérivées partielles suivante :

$$ \tilde{\mathcal{H}}_x u_x(y) = \bar{F}(x, y) - \bar{F}(x, y_x^*) \quad (2.33) $$

On suppose que l'équation (2.33) admet une solution $u_x \in C^1((\mathbb{R}^+)^d)$ vérifiant les propriétés suivantes :

* il existe une fonction $D_1 \in C^1((\mathbb{R}^+)^d)$ telle que :

$$ \forall \varepsilon > 0, \quad u_{x+\varepsilon} - u_x = \varepsilon D_1 $$

* il existe $p_0 \in \mathcal{P}$ et $c > 0$ tels que

$$ \forall y, \quad |\mathcal{H}_x u_x(x, y, \cdot)| \le c(V_{p_0} + 1) $$

où l'opérateur $\mathcal{H}_x$ est défini par (déjà défini par (2.9))

$$ \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) $$

Calculons $A^\varepsilon u_x(x, y, z)$ :

$$ \begin{aligned} A^\varepsilon u_x(x, y, z) = & \frac{r_0(y, z)}{\varepsilon} x [u_{x+\varepsilon}(y) - u_x(y)] + \frac{\mu_0(y, z)}{\varepsilon} x [u_{x-\varepsilon}(y) - u_x(y)] \\ & + \frac{\lambda_0(y, z)}{\varepsilon} x [u_{x-\varepsilon}(y + \varepsilon^p e_1) - u_x(y)] + \sum_{i=1}^{d} \frac{\mu_i(y, z)}{\varepsilon} y_i [u_x(y - \varepsilon^p e_i) - u_x(y)] \\ & + \sum_{i=1}^{d-1} \frac{\lambda_i(y, z)}{\varepsilon} y_i [u_x(y + \varepsilon^p (e_{i+1} - e_i)) - u_x(y)] \\ & + \frac{\lambda_d(y, z)}{\varepsilon} y_d [u_x(y - \varepsilon^p e_d) - u_x(y)] \end{aligned} $$

Comme $u_x$ est dérivable, en appliquant la formule de Taylor-Lagrange :

$$ u_x(y - \varepsilon^p e_i) - u_x(y) = -\varepsilon^p \frac{\partial u_x}{\partial y_i}(y) + R^1(y) $$
avec $R^1$ négligeable devant $\varepsilon^p$. En faisant de même avec tous les termes de $A^\varepsilon u_x(x, y, z)$ :

$$
\begin{align*}
A^\varepsilon u_x(x, y, z) ={}& r_0(y, z)xD_1(y) - \mu_0(y, z)xD_1(y) - \lambda_0(y, z)xD_1(y + \varepsilon^p e_1) \\
& + \frac{\lambda_0(y, z)}{\varepsilon^{1-p}}x\frac{\partial u_x}{\partial y_1}(y) - \sum_{i=1}^d \frac{\mu_i(y, z)}{\varepsilon^{1-p}}y_i\frac{\partial u_x}{\partial y_i}(y) \\
& + \sum_{i=1}^{d-1} \frac{\lambda_i(y, z)}{\varepsilon^{1-p}}y_i\left[\frac{\partial u_x}{\partial y_{i+1}}(y) - \frac{\partial u_x}{\partial y_i}(y)\right] - \frac{\lambda_d(y, z)}{\varepsilon^{1-p}}y_d\frac{\partial u_x}{\partial y_d}(y) \\
& + \frac{R^{\text{tot}}(y)}{\varepsilon}
\end{align*}
$$

avec le reste $R^{\text{tot}}$ négligeable devant $\varepsilon^p$.

Alors en multipliant $A^\epsilon u_x$ par $\epsilon^{1-p}$,

$$
\epsilon^{1-p} \times A^{\epsilon} u_{x}(x, y, z) = \mathcal{H}_{x} u_{x}(y) + o(1)
$$

Il nous reste donc à compenser l’erreur entre $\mathcal{H}_x u_x$ et $\tilde{\mathcal{H}}_x u_x$. Avec la même méthode que précédemment, comme cette erreur provient de la différence entre $Z^\epsilon$ et $\pi_{Y_d^\epsilon}$, nous introduisons à $x$ fixé l’équation de Poisson suivante associée à $(L_y, \mathcal{H}_x u_x)$ :

$$
\begin{align}
L_y h(z) &= \langle \mathcal{H}_x u_x(y, \cdot), \pi_y \rangle - \mathcal{H}_x u_x(y, z) \tag{2.34} \\
&= \tilde{\mathcal{H}}_x u_x(x, y) - \mathcal{H}_x u_x(x, y, z) \nonumber
\end{align}
$$

Avec la Proposition 6.2 et l'hypothèse 6.1, l'équation (2.34) ci-dessus a une solution $v_{x,y}$ bornée telle que pour tout $\epsilon > 0$, pour tout compact $K \in \mathbb{F}$ :

$$
\bullet \quad \forall (x, y, z) \in K, v_{x+\epsilon, y}(z) - v_{x,y}(z) = \mathcal{O}(\epsilon),
$$

$$
\bullet \quad \forall (x, y, z) \in K, v_{x, y+\epsilon}(z) - v_{x,y}(z) = \mathcal{O}(\epsilon).
$$

En choisissant $H_2^\epsilon(x, y, z) = \epsilon^{1-p}u_x(y) + \epsilon v_{x,y}(z)$, nous avons le résultat suivant :

**Proposition 6.4**

Pour tout compact $K \subset \mathbb{F}$,

$$
\lim_{\epsilon \to 0} \sup_{(x,y,z) \in K} \frac{1}{\sqrt{\epsilon}} \left| \bar{F}(x,y) - \bar{F}(x,y^*) - A^\epsilon H_2^\epsilon(x,y,z) \right| = 0
$$

**Démonstration :** En reprenant le développement de Taylor effectué sur le générateur ci-dessus, nous avons :

$$
\begin{align*}
A^{\epsilon} H_{2}^{\epsilon}(x, y, z) &= \mathcal{H}_{x} u_{x}(x, y, z) + L_{y} v_{x,y}(z) + \mathcal{O}(\epsilon^{p}) \\
&= \tilde{\mathcal{H}}_{x} u_{x}(x, y) + \mathcal{O}(\epsilon^{p}) \\
&= \bar{F}(x, y) - \bar{F}(x, y_{x}^{*}) + \mathcal{O}(\epsilon^{p})
\end{align*}
$$

car la fonction $u_x$ vérifie (2.33).

Alors on obtient le résultat car $p > \frac{1}{2}$ et que l'on travaille sur un compact $K$.

<figure><img src="image_28.png" /></figure>
### 6.2.4 Variation quadratique

La dernière étape consiste à calculer la variation quadratique du processus $X^\epsilon - H_{tot}^\epsilon$ afin de caractériser à l'aide du théorème central limite pour les martingales, la limite du terme (2.31).

**Proposition 6.5**

Soit pour tout $(x, y, z) \in \mathbb{F}$,

$$ 
\begin{aligned}
G(x, y, z) = x[r_0(y, z) + \mu_0(y, z)] &+ \left(\frac{\partial u_x}{\partial y_1}(y) - 1\right)^2 \lambda_0(y, z)x \\
&+ \sum_{i=1}^{d} \left(\frac{\partial u_x}{\partial y_i}(y)\right)^2 \mu_i(y, z)y_i + \sum_{i=1}^{d-1} \left(\frac{\partial u_x}{\partial y_{i+1}}(y) - \frac{\partial u_x}{\partial y_i}(y)\right)^2 \lambda_i(y, z)y_i \\
&+ \left[-(g_y(z + e_1) - g_y(z))x + (v_{x,y}(z + e_1) - v_{x,y}(z)) - \frac{\partial u_x}{\partial y_d}(y)\right]^2 \lambda_d(z)y_d \\
&+ \sum_{j=1}^{l} \left[(g_y(z - e_i) - g_y(z))x - (v_{x,y}(z - e_i) - v_{x,y}(z))\right]^2 \mu_{j+d}(z)z_j \\
&+ \sum_{j=1}^{l-1} \left[(g_y(z + e_{i+1} - e_i) - g_y(z))x - (v_{x,y}(z + e_{i+1} - e_i) - v_{x,y}(z))\right]^2 \lambda_{j+d}(z)z_j
\end{aligned}
$$

Alors,

$$ \frac{1}{\epsilon} [X^{\epsilon} - H_{tot}^{\epsilon}(\Psi^{\epsilon})]_{t} - \int_{0}^{t} G(\Psi^{\epsilon}) \, ds \xrightarrow{\epsilon \to 0} 0 $$

**Démonstration:** Pour illustrer le calcul, commençons par calculer le terme de $\frac{1}{\epsilon}[X^\epsilon - H_{tot}^\epsilon(\Psi^\epsilon)]_t$ lié au renouvellement $(x)$ devient $x + \epsilon$ avec un taux $xr_0/\epsilon$. En notant $P^{r_0}$ ce terme et avec $(x, y, z) = (X^\epsilon, Y^\epsilon, Z^\epsilon)$,

$$ P^{r_0}(x, y, z) = \int_{0}^{t} \left[ \epsilon - \epsilon^2 g_y(z) + \epsilon(v_{x+\epsilon,y}(z) - v_{x,y}(z)) + \epsilon^{1-p}(u_{x+\epsilon}(y) - u_x(y)) \right]^2 \frac{r_0(y, z)}{\epsilon^2} x \, ds $$

Comme $G$ représente le terme d'ordre 1 de la variation quadratique divisée par $\epsilon$, seuls les termes d'ordre $\epsilon$ à l'intérieur de la parenthèse nous intéressent. De plus, avec les hypothèses faites sur les dépendances des fonctions $u_x$, $v_{x,y}$ et $g_y$ (Hypothèses 6.1 et 6.2),

$$ P^{r_0}(x, y, z) = \int_{0}^{t} r_{0}(y, z) x \, ds + o(1) $$

Comme deuxième illustration, calculons le terme $P^{\mu_i}$ pour $1 \le i \le d$ :

$$ P^{\mu_i}(x, y, z) = \int_0^t \left[ -\epsilon(g_{y-\epsilon^p e_i}(z) - g_y(z))x + \epsilon(v_{x,y-\epsilon^p e_i}(z) - v_{x,y}(z)) \right. \\ \left. + \epsilon^{1-p}(u_x(y - \epsilon^p e_i) - u_x(y)) \right]^2 \frac{\mu_i(y, z)}{\epsilon^2} y_i \, ds $$

Au vu des dépendances des fonctions $g_y$ et $v_{x,y}$, le terme en violet est d'ordre $\epsilon^{1+p}$; et le terme restant s'écrit :

$$ u_x(y - \epsilon^p e_i) - u_x(y) = -\epsilon^p \frac{\partial u_x}{\partial y_i}(y) + o(\epsilon^p) $$
ce qui implique

$$P^{\mu_i}(x, y, z) = \int_0^t \left( \frac{\partial u_x}{\partial y_i}(y) \right)^2 \mu_i(y, z) y_i \, ds + o(1)$$

En poursuivant les calculs pour chaque transition possible du processus $\Psi^\epsilon$, on identifie la fonction $G$ voulue. ■

### 6.2.5 Passage à la limite

En assemblant tous les éléments des sections précédentes, nous en déduisons le théorème suivant :

#### Théorème 6.3

Pour tout $t > 0$, $M^\epsilon(t) = \frac{1}{\sqrt{\epsilon}}(M^{\epsilon,1}(t) - M^{\epsilon,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), y_{x(s)}^*, \cdot), \pi_{(y_{x(s)}^*)_d} \rangle \, ds$$

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

$$U^{\epsilon}(0) \xrightarrow{a.s} U^{in}$$

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

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

avec $W$ le mouvement brownien standard.

**Démonstration :** Soit $M^\epsilon$ la martingale suivante :

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

$H_{tot}^\epsilon$ est composé d'un terme d'ordre $\epsilon$ et d'un terme d'ordre $\epsilon^{1-p}$, ainsi

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

De plus avec la **Proposition 6.5** et d'après le théorème central limite appliqué à une martingale (Théorème 7.1.4 de [21]),

$$\begin{aligned} \mathbb{E}(M(t)M^T(t)) &= \lim_{\epsilon \to 0} \int_0^t G(X^\epsilon(s), Y^\epsilon(s), Z^\epsilon(s)) \, ds \\ &= \int_0^t \langle G(x(s), y_{x(s)}^*, \cdot), \pi_{(y_{x(s)}^*)_d} \rangle \, ds \end{aligned}$$

Enfin, en reprenant la décomposition de $U^\epsilon$ effectuée dans la **sous-section 6.2.1** et avec les résultats des sections précédentes, le processus $U^\epsilon$ converges vers le processus $U$ solution de

$$U(t) = U^{in} + \int_{0}^{t} [D'_{0}(x(s))x(s) + D_{0}(x(s))] U(s) \, ds + \int_{0}^{t} \sqrt{\langle G(x(s), y_{x(s)}^{*}, \cdot), \pi_{(y_{x(s)}^{*})_{d}} \rangle} \, dW_{s} \quad \blacksquare$$
## 6.3 Discussion de la limite obtenue

Le résultat obtenu dans le cas des trois échelles de temps s'interprète de la façon suivante :

* La partie déterministe du processus limite $U$ (Théroème 6.3) correspond à la dérivée de la fonction définissant l'EDO de la réserve $x$ du modèle limite. Ce terme était attendu puisque l'erreur est le terme suivant la limite dans le développement limité de $X^\epsilon : X^\epsilon = x + \sqrt{\epsilon U^\epsilon}$.

* La partie stochastique de $U$ est un mouvement brownien dont la variance est l'intégrale d'une fonction $G$ contre la limite $\Gamma$ de la mesure d'occupation de $(Y^\epsilon, Z^\epsilon)$ (voir le théorème 4.1). Chaque transition de $\Psi^\epsilon$ intervient dans cette variance. Cela est cohérent avec le fait que chaque transition possible ajoute de l'aléa au processus ce qui augmente la variance de celui-ci.

De plus, pour chaque transition agissant sur la réserve $X^\epsilon$ ($r_0$, $\mu_0$ et $\lambda_0$), la variance augmente du taux de la réaction correspondante. Les transitions agissant sur la partie intermédiaire $Y^\epsilon$ (($\mu_i$)$_{1\le i\le d}$ et ($\lambda_i$)$_{0\le i\le d}$) interviennent sur la variance via les dérivées premières de la fonction $u_x$. Cela est cohérent avec le fait que cette fonction est solution d'une équation sur le générateur $\tilde{H}_x$ (Hypothèse 6.2), qui est le générateur définissant la limite de la composante intermédiaire.

Enfin les transitions restantes (($\mu_{d+j}$)$_{1\le j\le l}$ et ($\lambda_{d+j}$)$_{0\le j\le l}$) concernent la partie rapide $Z^\epsilon$. Chacune de ces transitions est à l'origine d'un terme faisant intervenir les fonctions $g_y$ et $v_{x,y}$. Avec la même logique que précédemment, ces deux fonctions sont solutions d'une équation de Poisson (2.32) reliée à la distribution stationnaire $\pi$.

Ainsi chaque composante du processus limite (c'est-à-dire chaque échelle de temps) contribue à la variance de $U$.

Pour arriver à expliciter l'expression de l'erreur entre le modèle et sa limite, nous avons dû faire quelques hypothèses supplémentaires par rapport à celles faites pour prouver l'existence et l'unicité d'un modèle limite (Section 2).

Tout d'abord, nous avons restreint la vitesse de l'échelle intermédiaire en supposant que $p > 0.5$. Sans cette hypothèse, la fonction $u_x$ n'est plus suffisante pour contrôler l'erreur entre le taux de décroissance de la réserve moyennisé selon $\pi$ et le même taux évalué au point fixe $y^*$ (terme (2.29)). L'adaptation de cette méthode au cas $p \le 0.5$ est un problème intéressant, nécessitant probablement l'introduction de nouveaux objets mathématiques voire une autre décomposition de l'erreur afin de pouvoir contrôler le terme issu de l'échelle intermédiaire.

Pour effectuer le calcul de la variation quadratique de la martingale $M^\epsilon$, il est nécessaire de connaître la dépendance en $x$ et/ou $y$ des fonctions solutions des différentes équations introduites (une équation aux dérivées partielles et deux équations de Poisson). Pour ce faire, en nous inspirant du résultat obtenu dans le cas du modèle à 2 échelles de temps, nous avons supposé une dépendance uniforme et linéaire (hypothèses 6.1 et 6.2). Une étude un peu plus précise des solutions de l'équation de Poisson, tirant parti du fait que le générateur $L_y$ est linéaire, permettrait probablement de prouver ces hypothèses. Ce travail est en cours.

Pour l'équation aux dérivées partielles (2.33), la situation est un peu plus compliquée. Nous n'avons pas réussi à prouver l'existence de solution malgré la simplicité de cette équation linéaire aux dérivées partielles du premier ordre. L'étude de la dépendance des solutions n'étant ainsi pas envisageable, nous avons fait la même hypothèse que pour l'équation de Poisson.
Cette étude de l’erreur se focalise sur la réserve. Comme celle-ci a été renormalisée par $\varepsilon$ pour tenir compte du très grand nombre de follicules au sein de la réserve, si l’on effectue un développement limité sur les termes du générateur dépendant de la réserve (Définition 2.1), le fait de considérer l’erreur à l’ordre $\sqrt{\varepsilon}$ devient évident. Ainsi la définition de $U^\varepsilon$ n’est pas anodine et fait partie intégrante de la preuve.

L’étude d’autres types d’erreur (sur les mesures d’occupation de $Y^\varepsilon$ ou de $Z^\varepsilon$) nécessiterait probablement de diviser la différence entre le modèle et sa limite par une autre puissance que $\sqrt{\varepsilon}$. Un changement d’ordre de l’erreur pourrait aussi être une piste prometteuse pour adapter la méthode utilisée au cas $p \le 0.5$.

Enfin d’un point de vue plus pratique, le modèle limite obtenu semble très peu utilisable de part sa complexité : il fait intervenir trois fonctions qui sont intégrées contre la distribution $\pi$. Malgré tout, au moins numériquement, le processus $U$ peut être simuler ce qui pourrait permettre de mieux comprendre les hypothèses implicites utilisées si l’on confronte le modèle limite à des données expérimentales.
<figure><img src="image_29.png" /></figure>
# Chapitre 3

## Comparison of the model with experimental data on ovarian follicle population

This last chapter presents an ODE-based model adapted to the biological experimental data available on follicles numbers in the mouse. Taking into account a study of the theoretical identifiability of the parameters, we design and perform an estimation of the parameters of the model. Then, we illustrate the model potential by reproducing in silico experiments.

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

## Sommaire

<table><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><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></table>
<table>
  <tbody>
    <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>
  </tbody>
</table>
# 1 Introduction

The ovary is a unique instance of dynamic, permanently remodeling organ in adulthood. Such a remodeling results from the asynchronous development of ovarian follicles, which starts with activation from a quiescent pool and ends up either by ovulation or degeneration. Ovarian follicles are somatic structures sheltering the oocyte (female germ cell) and secreting a variety of hormones and growth factors; they underlie both the endocrine and exocrine (gamete release) facets of ovarian function. The physiological status of the ovaries is determined by the distribution of follicles into different maturity stages and by the total number of follicles [11]. In mammals, females are endowed with a pool of quiescent (primordial) follicles established early in life, whose progressive exhaustion is the prominent marker of ovarian aging [13, 8].

Despite the critical importance of assessing accurately the number of follicles, very few datasets are available to follow the follicle populations on a whole life scale. In vivo, only the terminal stages can be monitored through ultrasonography imaging in large species [1, 3]. To detect all follicles, invasive approaches, relying on histology are needed. Ovaries are fixed, serially sectioned and stained with proper dyes, and manually analyzed by light microscopy. Such a counting is a tedious and time consuming procedure requiring cutting-edge expertise. As a consequence, although the great majority of studies dedicated to female reproductive biology provides one with partial information of follicle numbers, only a couple of them is suitable for modeling follicle population dynamics and estimating the growth and death rates affecting the follicles according to the development stage.

A series of studies from J.M. Faddy and co-authors have been dedicated to fitting compartmental linear models on datasets providing the number of follicles according to age, especially in the mouse [24, 23, 25]. To cope with the changes in follicle populations, they needed to introduce piecewise constant transition rates and choose empirically the switching time values. At the time these studies were performed, age-varying rates were a natural choice since very little was known on the control of follicle activation and early growth. Since then, biological knowledge has been greatly improved, especially thanks to the use of mouse genetic models. It is now clear that growing follicles affect the activation rate through secreting the anti-Müllerian hormone (AMH) that acts locally on quiescent follicles [19, 18]. From the modeling viewpoint, this means that interactions between follicles cannot be ignored and that a nonlinear framework is more appropriate for models of follicle populations.

In addition to knock-out (KO) models, the use of transgenic mice to label and image specifically somatic and/or germ cells, in a time-controlled manner, have shed new light onto the formation of ovarian follicles. Besides revealing antero-posterior waves and a dorso-ventral gradient during ovarian development [12, 62], these studies have confirmed pioneering observations on the heterogeneity of the follicle pool in early life [32]. It appears that two distinct populations of follicles coexist early in life. One emanates from the inner part of the ovaries (the medulla) before birth, while the other is formed during the first postnatal days in the outer part (cortex). The latter will build up the ovarian reserve and sustain the reproductive function during adulthood. The former, usually designated as the first wave of follicles, enters growth almost directly, participates in the very first ovarian cycles and vanishes in a few weeks [50, 66]. From a modeling viewpoint, the shift in the onset of the first wave compared to the reserve, coupled with the timescale difference in the activation rate is susceptible to affect significantly the dynamics of the whole follicle population according to age in the first weeks.

In this work, we first design an ODE-based data-driven model, in which the number of compartments is dictated by the developmental stages considered in the available datasets on follicle numbers according to age. Following the framework introduced in [6] and analyzed in the stochastic
tic case in [4], we account for the interactions between follicle stages by means of a nonlinear activation rate. We then complete this core model to investigate more deeply the transient period when both the first wave and the reserve of follicles coexist. We add one compartment upstream from the less mature stage to capture the influence of the rate of follicle formation and thereafter duplicate the number of compartments to follow separately the dynamics of the follicle pool emanating from the first wave. We perform a theoretical identifiability analysis, and then settle a parameter estimation strategy based on a cost function merging the different sources of experimental data: total follicle numbers in each follicle stage, information from AMH-KO mice in which the activation rate is unregulated, and information on the relative proportion of the first wave follicles among the growing follicle pool. The results highlight the interest of accounting for prefollicular steps in oocyte development to monitor accurately the follicle dynamics on the long-term, and possibly forecast the status of the follicle reserve with ovarian aging. They also give cues to investigate separately the growth and early extinction of the first wave-derived follicle pool, which are still poorly understood despite their considerable importance in programming adult reproductive life.

The paper is organized as follows. In the second section, we describe the different datasets that will be used for the model design and parameter estimation. In the third section, we introduce the core model and the two-pool model, and analyze the theoretical identifiability of the model. The fourth section details the optimization strategy and statistical regression procedure, including statistical comparisons of nested models. We finally exploit the estimation results to suggest targeted additional information and illustrate predictions of *in silico* experiments (fifth section), before ending by a general discussion section.

# 2 Experimental datasets on ovarian follicle populations

## 2.1 Longitudinal data on ovarian follicles numbers as a function of age

In this study, we use the datasets on mouse ovarian follicles provided in [23, 25]. These datasets distribute the total number of germ cells reported in the original study of [34] (see the bottom right panel of Figure 3.1) into five different classes inspired from the classification introduced in [45]. These classes are based on morphological criteria affecting the somatic cells surrounding the oocyte, namely the cell shape, number of cell layers and possible presence of a cavity (antrum) (see Table 3.1).

<table><caption>TABLE 3.1 – Classification used in the datasets available in [23, 25].</caption><thead><tr><td>Class</td><td>Description</td></tr></thead><tbody><tr><td>I</td><td>primordial follicles: single layer of squamous pregranulosa cells</td></tr><tr><td>II</td><td>primary follicles: single layer of cuboidal granulosa cells and growing oocyte</td></tr><tr><td>III</td><td>secondary follicles with two cell layers</td></tr><tr><td>IV</td><td>secondary follicles with three cell layers</td></tr><tr><td>V+</td><td>all other growing follicles: four or more cell layers, with or without an antrum</td></tr></tbody></table>

It is worth noting right now that the classification departs from the original proposition in [45], since antral follicles are pooled with pre-antral ones in a single class V+, which will impact the formulation of the model. Compared to the widely-used classification established in the seminal work of Pedersen [56], an additional difference lies in the classification of the most immature stages. Again, Class I merges two types of follicles, the quiescent (primordial) follicles and transitional
follicles, which are undergoing activation. The separation between both classes is subtle from a morphological viewpoint, yet possibly important for a study focused on follicle activation, as the merging might overassess the number of fully quiescent follicles.

Dataset A₁ [25] provides the number of ovarian follicles from birth (namely from one day after birth) until an age of 100 days, with an age-resolution of one day in a single mouse strain (CBA/A).

Dataset A₂ [23] provides the number of ovarian follicles on a whole life range (about 500 days) with an age-resolution of about 8 days in 4 mouse strains, including the CBA strain. Thanks to the genetic proximity of the CBA strain with the CBA/A strain used in [25], we can combine dataset A₁ and dataset A₂ in a single dataset (Dataset A) (see respectively the blue and orange datapoints on Figure 3.1).

<figure><img src="image_30.png" /><figcaption>**FIGURE 3.1** – Follicle numbers as a function of age (Dataset A).
Dataset A combines Dataset A₁ (blue points) extending from birth to 100 days with Dataset A₂ (orange points) covering the whole life. From left to right and top to bottom : follicle numbers in class I, II, III, IV and V. The rightmost bottom panel shows the total number of follicles (purple points) in all 5 classes.</figcaption></figure>

**2.2 Additional specifications**

**Control of the activation rate** As stated in the Introduction, the contribution of AMH-mediated interactions between follicles on the activation rate can be assessed from AMH-KO mice. The comparative study performed in [19] provides us with additional data point at ages 25 and 120 days, on the numbers of both quiescent (primordial) follicles and growing follicles in either control mice or homozygous AMH-/- mice.

Although these data were acquired in strain C57BL/6 (the most frequently used strain in current reproductive biology), the control datapoints fall within the variability range of Dataset A. We can thus take advantage of this study (Dataset B) to compare the control and perturbed situation which results in a accelerated depletion of the quiescent follicle pool (compare the black bars to the blue ones on the left panel of Figure 3.2).

**Two populations of follicles** The contribution of the two distinct follicle populations co-existing early in life has been quantified in [66], through the proportion of first-wave derived follicles among all growing follicles (right panel of Figure 3.2). Note that there is no standard deviation for the first and the last three points.
<figure><img src="image_31.png" /><figcaption>FIGURE 3.2 – Additional experimental data.</figcaption></figure>

Left panel. Dataset *B*: Comparison of follicle numbers in wild type mice (black bars) and AMH-KO mice (blue bars) at two different ages (left column: 25 day old mice, right column: 120 day old mice).

Right panel. Dataset *C*: Proportion (in percent) of growing follicles of the first wave among all growing follicles. For both Datasets *B* and *C*, the vertical red bars represent the standard deviation.

No digital data were available; we proceeded manually by clicking the datapoint coordinates directly on the figures for Datasets *A* and *B*, and retrieved raw numbers from the tables in [66] for Dataset *C*.

# 3 Compartmental model of ovarian follicle dynamics with aging

## 3.1 Model design

The model consists of two layers.

The core layer is data-driven; the number of compartments and their interpretation as follicle maturation stages is directly dictated by the available datasets, and primarily by Dataset *A*. The core model consists of 5 compartments, one dedicated to the quiescent follicles (compartment I) and the 4 other ones (compartments II to V) to the successive growing stages as described in Table 3.1. In each compartment, ovarian follicles are subject to two types of event; they can either move forward to the next compartment (unidirectional migration, except in the last compartment) or degenerate and disappear (death). The activation rate from the first compartment is controlled by the compartments corresponding to the secondary follicle stage, which contributes the most to AMH-secretion [51], namely compartments III and IV.

To capture more accurately the follicle dynamics during the first days of life, the core model is embedded into a more detailed model (see Figure 3.3). First, to deal with the early step of follicle formation, we add a compartment (compartment 0) corresponding to the pre-follicular stage of oogenesis (where oocytes are gathered in syncytium structures, the germ cell cysts). Germ cells may self-renew, undergo apoptosis, or move to compartment I. The transition from the pre-follicular compartment to the next one corresponds to the formation of ovarian follicles. Then, to account for
the existence of two distinct populations of follicles, the model is duplicated. We end up with two population models running almost in parallel, with population-specific rates of migration and birth. Coupling between the dynamics occurs at the level of the activation rate of quiescent follicles from the reserve, since all AMH-secreting follicles participate in the feedback exerted on this rate. Follicles from the first wave begin to grow as soon as they are formed. Accordingly, the transition rate between compartments I and II is not controlled in the model, and the whole dynamics of the first wave is linear.

<figure><img src="image_32.png" /><figcaption>FIGURE 3.3 – Schematic view of the compartmental model (S).</figcaption></figure>

The orange box represents the core 5-compartment model. The green box represents the pre-follicular stage. The blue box represents the 6 compartments structuring the first wave dynamics.

Let, for $0 \le i \le 5$, $\bar{y}_i(t)$ and $y_i(t)$ be the number of follicles in compartment $i$ at time $t$, and emanating respectively from the first wave or the reserve. We formulate the model as the following system of ordinary differential equations (using the same color code as in Figure 3.3):

$$ (S): \left\{ \begin{aligned} \frac{d\bar{y}_0}{dt}(t) &= [\bar{r}_0 - (\bar{m}_0 + \bar{g}_0)]\bar{y}_0 \\ \frac{d\bar{y}_i}{dt}(t) &= \bar{m}_{i-1}\bar{y}_{i-1} - (\bar{m}_i + \bar{g}_i)\bar{y}_i, \quad \text{for } 1 \le i \le 5 \\ \frac{dy_0}{dt}(t) &= [r_0 - (m_0 + g_0)]y_0 \\ \frac{dy_1}{dt}(t) &= m_0y_0 - [\lambda_1(\bar{y}(t), y(t), t) + g_1]y_1 \\ \frac{dy_2}{dt}(t) &= \lambda_1(\bar{y}(t), y(t), t)y_1 - (m_2 + g_2)y_2 \\ \frac{dy_i}{dt}(t) &= m_{i-1}y_{i-1} - (m_i + g_i)y_i, \quad \text{for } 3 \le i \le 5 \end{aligned} \right. \quad (3.1) $$

with $\bar{m}_5 = m_5 = 0$ and $\lambda_1(\bar{y}, y, t)$ the controlled activation rate of quiescent follicles:

$$ \lambda_1(\bar{y}, y, t) = \underbrace{\left[ m_1 + \frac{f_1}{1 + [K_1(\bar{y}_3 + y_3 + \bar{y}_4 + y_4)]^n} \right]}_{\text{control term}} \times \underbrace{\frac{(t/\tau)^p}{1 + (t/\tau)^p}}_{\text{delay term}} \quad (3.2) $$

The activation rate of quiescent follicles combines a control term, representing the feedback of the growing, AMH-secreting follicles, with a delay term, accounting for the shift observed in early life in the onset of follicle formation and activation in the reserve follicle pool compared to the
first-wave pool. For $p$ large enough (we set for instance $p = 20$), the delay function behaves as the unity step function : it is closed to zero for small $t$ ($t < \tau$) and closed to 1 when $t$ exceeds $\tau$. All rates are strictly positive except the renewal and death rates ($\bar{r}_0, r_0, \bar{g}_i$ and $g_i$), which may be zero. The initial conditions are the following :

$$ \left\{ \begin{array}{l} \bar{y}_0(t_0) = \bar{y}_{init} > 0, \\ y_0(t_0) = y_{init} > 0, \\ \bar{y}_i(t_0) = y_i(t_0) = 0, \text{ for } i \ge 1. \end{array} \right. $$

To represent the dynamics of the follicle population right from the first day after birth (earliest age in Dataset A), we have to choose an initial time during embryonic life. We set $t_0 = E16.5$, hence 3 days before birth (which occurs at embryonic age E19.5), just after the peak reached by the germ cells [52] (from this time on the germ cell number decreases) and close to the appearance of the very first follicles in the dorsal medulla area of the ovaries [12]. Since time represents postnatal age in the model, note that $t_0 = -3d$. During the period spanning a few days before birth (E16.5-E19.5) until one week of (postnatal) age, the population of germ cells is subject to a series of events happening sequentially on the individual level but partly overlapping, due to between-cell desynchronization, on the population level : ending of proliferative rounds, cell cycle exit, entry into meiosis more or less concomitant with embedding within follicles, and associated with a high risk of cell death [54]. These events are represented in a compact way in the germ cell dynamics $y_0(t)$ and $\bar{y}_0(t)$, which includes together self-renewal at rate $r_0/\bar{r}_0$, cell death at rate $g_0/\bar{g}_0$ and formation of follicles at rate $m_0/\bar{m}_0$. Desynchronization exists both between the subpopulations (with a clear advance of the first wave) and within both subpopulations (for instance the formation of the follicle reserve is spread over several days).

Note that the state variables of the complete model are augmented with respect to the core layer, without significantly increasing the number of observable states. We can only observe the sum of the follicle numbers ($y_i^{tot}(t) = \bar{y}_i(t) + y_i(t)$) from both populations in compartments I to V from Dataset A, and get information on the proportion of growing follicles emanating from the first wave from Dataset C. No direct observation is usable for the pre-follicular compartments, which will ground the identifiability analysis performed next in Sous-section 3.2.

## 3.2 Model analysis : theoretical identifiability

In this section, we first reparameterize the model, define the observable variables and introduce the notions of local and global identifiability. Then, we study separately the identifiability of the model parameters and initial conditions, after introducing a relevant reduced model. Finally we conclude on the identifiability of the whole system.

### 3.2.1 Model reparameterization and reduced model

We first introduce the notations used for this study.

$$ \text{Let} \quad \left\{ \begin{array}{l} \bar{D}_0 = \bar{g}_0 + \bar{m}_0 - \bar{r}_0 \\ \bar{D}_i = \bar{m}_i + \bar{g}_i, \text{ for } 1 \le i \le 5 \\ D_0 = m_0 + g_0 - r_0 \\ D_1 = g_1 \\ D_i = g_i + m_i, \text{ for } 2 \le i \le 5 \end{array} \right. \quad \text{and} \quad \left\{ \begin{array}{l} \bar{Y}_0(t) = \bar{m}_0 \bar{y}_0(t) \\ Y_0(t) = m_0 y_0(t) \end{array} \right. , $$

then model (S) ((3.1)) becomes :
$$ (S1): \left\{ \begin{array}{l} \frac{d\bar{Y}_0}{dt}(t) = -\bar{D}_0\bar{Y}_0 \\ \frac{d\bar{y}_1}{dt}(t) = \bar{Y}_0 - \bar{D}_1\bar{y}_1 \\ \frac{d\bar{y}_i}{dt}(t) = \bar{m}_{i-1}\bar{y}_{i-1} - \bar{D}_i\bar{y}_i, \text{ for } 2 \le i \le 5 \\ \frac{d\bar{Y}_0}{dt}(t) = -D_0Y_0 \\ \frac{d\bar{y}_1}{dt}(t) = Y_0 - D_1y_1 - \lambda_1(y, \bar{y}, t)y_1 \\ \frac{d\bar{y}_2}{dt}(t) = \lambda_1(y, \bar{y}, t)y_1 - D_2y_2 \\ \frac{d\bar{y}_i}{dt}(t) = m_{i-1}y_{i-1} - D_iy_i, \text{ for } 3 \le i \le 5 \end{array} \right. , $$

$$ \text{with } \lambda_1(y, \bar{y}, t) = \left[ m_1 + \frac{f_1}{1 + [K_1(y_3^{\text{tot}}(t) + y_4^{\text{tot}}(t))]^n} \right] \left( \frac{(t/\tau)^{20}}{1 + (t/\tau)^{20}} \right). $$

The new parameter vector $\theta^I$ reads :

$$ \theta^I = \left[ (\bar{D}_i)_{0 \le i \le 5}, (\bar{m}_j)_{1 \le j \le 4}, (D_i)_{0 \le i \le 5}, (m_j)_{1 \le j \le 4}, f_1, K_1, n, \tau, \bar{Y}_0(t_0), (\bar{y}_i(t_0))_{1 \le i \le 5}, Y_0(t_0), (y_i(t_0))_{1 \le i \le 5} \right] \quad (3.3) $$

With this new parameterization, $\bar{m}_0$ and $m_0$ are embedded into the initial conditions (resp. $\bar{Y}_0(t_0), Y_0(t_0)$).

Given the available datasets described in Section 2 and the different submodels involved in this study, we define the following observable variables for any $t \ge t_0$,

$$ \left\{ \begin{array}{l} O_i^{\text{tot}}(t) = \bar{O}_i(t) + O_i(t) = y_i^{\text{tot}}(t), \quad \text{for } 1 \le i \le 5, \quad \text{(Dataset A)} \\ P_0(t) = 100 \times \frac{\sum_{i=2}^{5} \bar{y}_i(t)}{\sum_{i=2}^{5} y_i^{\text{tot}}(t)}, \quad \text{(Dataset C)} \\ \bar{O}_i(t) = \bar{y}_i(t), \quad \text{for } 1 \le i \le 5, \\ O_i(t) = y_i(t), \quad \text{for } 1 \le i \le 5, \end{array} \right. $$

Note that we do not consider Dataset B for the theoretical identifiability, since the AMH-KO condition corresponds to a simpler model with uncontrolled activation rate $\lambda_1$ ($K_1 = 0$).

Also, the observables are supposed to be available at any time $t \ge t_0$, while the experimental observations are only available at point time values, and on a finite time horizon.

Any parameter $\theta_j^I$ is said to be

* globally identifiable, if every solution of (S1) consistent with the set of observables ($P_0, (O_i^{\text{tot}})_{1 \le i \le 5}$) involves the same value $\theta_j^I$,

* locally identifiable, if there is a finite number of $\theta_j^I$ values among the solutions of (S1) consistent with the set of observables ($P_0, (O_i^{\text{tot}})_{1 \le i \le 5}$),

* not identifiable, if there is an infinite number of $\theta_j^I$ values among the solutions of (S1) consistent with the set of observables ($P_0, (O_i^{\text{tot}})_{1 \le i \le 5}$).

**Reduced model** In order to reduce the number of states and their related parameters, we start with removing the last compartments, whose dynamics do not impact the model outputs in the previous compartments. We thus analyze subsystem (S2) with state variables : ($\bar{Y}_0, \bar{y}_1, \bar{y}_2, \bar{y}_3, \bar{y}_4, Y_0, y_1, y_2, y_3, y_4$) with the corresponding subset of observables ($O_i^{\text{tot}})_{1 \le i \le 4}$.
### 3.2.2 Identifiability of the model parameters of the reduced model

In a first step, we only analyze the identifiability of the kinetic parameters (initial conditions are left out), using the *Structural Identifiability* Julia package [16]. This package cannot handle non-autonomous systems, neither rate functions with a non-integer power. To overcome these limitations, and get information on the delay $\tau$ and exponent $n$ appearing in (3.2), we introduce:

* $\frac{dw_1(t)}{dt} = 1$, or, equivalently, $w_1(t) = (t - t_0) + w_1(t_0)$,

* $\frac{dw_2(t)}{dt} = n \frac{w_2}{y_3^{tot} + y_4^{tot}} (y_3^{tot} + y_4^{tot})'$,

or, equivalently, $w_2(t) = [K_1(y_3^{tot}(t) + y_4^{tot}(t))]^n + C$ with $C = w_2(t_0) - [K_1(y_3^{tot}(t_0) + y_4^{tot}(t_0))]^n$.

Note that $K_1$ is embedded into the initial condition $w_2(t_0)$.

Then, we have $\lambda_1(y, \bar{y}, t) = \tilde{\lambda}_1(w_1, w_2) = [m_1 + \frac{f_1}{1 + w_2}] (\frac{(w_1/\tau)^{20}}{1 + (w_1/\tau)^{20}})$, and we can write the autonomous system (S2') equivalent to (S2), with state variables

$$ (S2') : \begin{cases} \bar{Y}_0, & (\bar{y}_i)_{1\le i \le 4} \\ Y_0, & (y_i)_{1\le i \le 4} \\ w_1, & w_2 \end{cases} \qquad (3.4) $$

Since $w_1$ represents time and is defined up to a constant, the observable variables are $(O_j^{tot})_{1\le j \le 4}$, and $w_1$.

Due to the number of states, the identifiability analysis is hardly computationally tractable. To save computational time, we use the following trick, which was suggested to us by Gleb Pogudin<fnref n="1" /> (personal communication).

The trick is based on augmenting the dimension of both the state space and observable space. We introduce the integrals of the observable variables $(O_j^{tot})_{1\le j \le 4}$:

$$ I'_{j}(t) = O_{j}^{tot}(t), \text{ for } 1 \le j \le 4 $$

and let:

$$ \left\{ \begin{aligned} J_1 &= I_2 + \frac{D_1}{D_2} I_1 + \frac{1}{D_2} \left(y_1 + y_2 + \frac{Y_0}{D_0}\right) + \frac{1}{D_2} \bar{y}_2 + \frac{\bar{D}_2 D_1 + \bar{m}_1 D_2}{\bar{D}_1 \bar{D}_2 D_2} \left(\bar{y}_1 + \frac{\bar{Y}_0}{\bar{D}_0}\right) \\ J_2 &= I_3 + \frac{D_1 m_2}{D_2 D_3} I_1 + \frac{m_2}{D_2 D_3} \left(\frac{Y_0}{D_0} + y_1 + y_2 + \frac{D_2}{m_2} y_3\right) + c_1 \left(\bar{y}_1 + \frac{\bar{Y}_0}{\bar{D}_0}\right) + \frac{\bar{m}_2}{D_2 D_3} \left(\bar{y}_2 + \frac{\bar{D}_2}{\bar{m}_2} \bar{y}_3\right) \\ J_3 &= I_4 + \frac{m_2 m_3}{D_2 D_3 D_4} \left[D_1 I_1 + \frac{Y_0}{D_0} + y_1 + y_2 + \frac{D_2}{m_2} y_3 + \frac{D_2 D_3}{m_2 m_3} y_4\right] + c_2 \left(\bar{y}_1 + \frac{\bar{Y}_0}{\bar{D}_0}\right) \\ &\quad + \frac{\bar{m}_3}{D_3 D_4} \left(\frac{\bar{m}_2}{\bar{D}_2} \bar{y}_2 + \bar{y}_3 + \frac{\bar{D}_3}{\bar{m}_3} \bar{y}_4\right) \end{aligned} \right. \qquad (3.5) $$

with $c_1 = \frac{D_1 m_2}{D_1 D_2 D_3} + \frac{\bar{m}_1 \bar{m}_2}{\bar{D}_1 \bar{D}_2 \bar{D}_3}$ and $c_2 = \frac{D_1 m_2 m_3}{D_1 D_2 D_3 D_4} + \frac{\bar{m}_1 \bar{m}_2 \bar{m}_3}{\bar{D}_1 \bar{D}_2 \bar{D}_3 \bar{D}_4}$.

We design linear combinations between the $I_i$'s and state variables (3.5), in such a way that $J_1, J_2$ and $J_3$ are constant. As a result, we can express $I_2, I_3$ and $I_4$ in terms of $I_1, J_1, J_2, J_3$, and replace states $I_2, I_3$ and $I_4$ by the new parameters $J_1, J_2$ and $J_3$.

<footnote n="1">1. LIX, CNRS, École Polytechnique, Institut Polytechnique de Paris, France gleb.pogudin@polytechnique.edu</footnote>
This first step leads to the following conclusions regarding the identifiability of the parameters of system (S2) :

• $(\bar{D}_i)_{0 \le i \le 4}, (\bar{m}_j)_{1 \le j \le 3}, (D_i)_{0 \le i \le 4}, (m_j)_{2 \le j \le 3}$ and $(J_i)_{1 \le i \le 3}$ are globally identifiable,

• $m_1, f_1, n$ and $\tau$ are locally identifiable.

### 3.2.3 Identifiability of the initial conditions of the reduced model

In a second step, we analyze the identifiability of the initial conditions of system (S2), using the SIAN Julia package 2 (also implemented in Maple 3) [33]. Using the results of the previous subsection, the globally identifiable parameters are added to the observables, and we obtain that:

• $(\bar{y}_i(t_0))_{1 \le i \le 4}, (y_i(t_0))_{1 \le i \le 4}, \bar{Y}_0(t_0), Y_0(t_0)$ and $w_1(t_0)$ are globally identifiable,

• $w_2(t_0)$ is locally identifiable.

### 3.2.4 Identifiability of the initial system

Given that the parameters ruling the dynamics of the follicles emanating from both the first wave and the reserve are identifiable, we can finally use separately the $\bar{O}_i(t)$ and $O_i(t)$ observables (instead of the aggregate observations $O_i^{tot}$), (for $1 \le i \le 4$), together with $P_0(t)$, to study the identifiability of system (S1), and conclude that:

• $m_4, m_5, D_4$ and $D_5$ are globally identifiable,

• the initial conditions $(\bar{y}_5(t_0)$ and $y_5(t_0))$ are globally identifiable.

Going back to the initial parameters of system (S) ((3.1)), we can sort the model parameters as follows:

• globally identifiable parameters : $(\bar{g}_i)_{1 \le i \le 5}, (\bar{m}_j)_{1 \le j \le 4}, (g_i)_{1 \le i \le 5}$ and $(m_j)_{2 \le j \le 4}$,

• locally identifiable parameters : $m_1, f_1, K_1, n$ and $\tau$,

• non identifiable parameters : $\bar{g}_0, \bar{m}_0, \bar{r}_0, g_0, m_0, r_0, \bar{y}_{init}$ and $y_{init}$.

One can note that the parameters related to the germ cell compartments are not identifiable, which is rather intuitive since there are no information on the germ cell compartments in the observable datasets. Also, we cannot distinguish the death rate from the renewal rate in compartment 0. Nevertheless, including the germ cell dynamics into the model improves the monitoring of follicle dynamics in early life, as detailed in Section 4. In addition, we can take advantage of the identifiability of the compound parameters to get insight into the germ cell dynamics : the timescale parameters $(\bar{D}_0 = \bar{g}_0 + \bar{m}_0 - \bar{r}_0$ and $D_0 = m_0 + g_0 - r_0)$ and the initial rate of quiescent follicle formation $(\bar{Y}_0(t_0) = \bar{m}_0\bar{y}_0(t_0)$ and $Y_0(t_0) = m_0y_0(t_0))$ can be identified, as will be developed in Sous-section 5.1.

<footnote n="2">2. https://github.com/alexeyovchinnikov/SIAN-Julia</footnote>

<footnote n="3">3. https://maple.cloud/app/6509768948056064/Structural+Identifiability+Toolbox</footnote>
# 4 Parameter estimation strategy

In this section, we design an optimization strategy suited for handling simultaneously the different datasets. We first adimension most of the parameters, and then define relevant model errors for each data type. We optimize the composite cost function and check the practical identifiability of the estimated parameters. Finally, we perform a statistical model selection and compare the full original model with simpler models to avoid data overfitting.

## 4.1 Parameter renormalization and search intervals

To deal with as many dimensionless parameters as possible, and define intrinsic bounds on the parameter values, we reparametrize the model as follows:

$$ \bullet \bar{l}_0 = \frac{\bar{g}_0 + \bar{m}_0 - \bar{r}_0}{g_0 + m_0 - r_0} $$

$$ \bullet \bar{l}_1 = \frac{\bar{m}_1 + \bar{g}_1}{l_1} $$

$$ \bullet \bar{l}_i = \frac{\bar{m}_i + \bar{g}_i}{m_i + g_i} \text{ for } 2 \le i \le 5, $$

$$ \bullet \bar{\rho}_i = \frac{\bar{m}_i}{\bar{m}_i + \bar{g}_i} \text{ for } 1 \le i \le 4, $$

$$ \bullet \tilde{f}_1 = \frac{f_1}{m_1 + f_1} $$

$$ \bullet l_0 = \frac{g_0 + m_0 - r_0}{l_1} $$

$$ \bullet l_i = \frac{m_i + g_i}{l_1} \text{ for } 2 \le i \le 5, $$

$$ \bullet \rho_1 = \frac{m_1 + f_1}{m_1 + f_1 + g_1}, $$

$$ \bullet \rho_i = \frac{m_i}{m_i + g_i} \text{ for } 2 \le i \le 4, $$

Parameter $l_1 := m_1 + f_1 + g_1$, the decay rate of the quiescent follicles in the reserve, is chosen as the reference timescale for all compartments. We only consider ratios of parameters, among which some range naturally between 0 and 1 (see Table 3.2). The remaining parameters entering the activation rate, $\tau$, $n$, and $K_1$ are kept unchanged.

In addition, we will estimate reparameterized initial conditions

$$ \bullet \hat{N}_{tot} = \frac{\bar{Y}_0(t_0) + Y_0(t_0)}{l_1} = \frac{\bar{m}_0 \bar{y}_0(t_0) + m_0 y_0(t_0)}{l_1}, \text{ the rescaled total number of germ cells at initial time,} $$

$$ \bullet \hat{r}_{init} = \frac{\bar{Y}_0(t_0)}{\bar{Y}_0(t_0) + Y_0(t_0)} = \frac{\bar{m}_0 \bar{y}_0(t_0)}{\bar{m}_0 \bar{y}_0(t_0) + m_0 y_0(t_0)}, \text{ the rescaled initial proportion of germ cells committed to the first wave.} $$

The new parameter vector $\theta$ reads:

$$ \theta = \left[ (\bar{l}_i)_{0 \le i \le 5}, (\bar{\rho}_j)_{1 \le j \le 4}, (l_i)_{0 \le i \le 5}, (\rho_j)_{1 \le j \le 4}, \tilde{f}_1, K_1, \tau, n, \hat{N}_{tot}, \hat{r}_{init} \right] \quad (3.6) $$

Table 3.2 gathers the bounds and scale (log or linear scale) used for parameter search. The bounds have been defined according to the following biological specifications:

1. The first wave follicles are only involved in the first ovarian cycles, contrary to follicles emanating from the reserve, which cover the whole reproductive life. The dynamics of the first wave is faster, so that for any $0 \le i \le 5$, $\bar{l}_i$ (the normalized exit rate from compartment i for the first wave) is lower bounded by 1.

2. For the reserve decay rate, $l_1$, we choose a search interval $[10^{-4}, 10^0]$ two logs above and under the order of magnitude of the half-life of the reserve (on the order of $10^2$, as estimated
<table><caption>**TABLE 3.2** – Estimation bounds and use of $log_{10}$ scale for the model parameters

Parameters above the double horizontal bar are searched in $log_{10}$ scale, those below in linear scale. The values fixed from biological specifications are highlighted in red, the intrinsic bounds in green.</caption><thead><tr><th></th><th>Bounds</th></tr></thead><tbody><tr><td>(l<sub>i</sub>)<sub>0&le;i&le;5</sub></td><td>[1, 10<sup>3</sup>]</td></tr><tr><td>l<sub>0</sub></td><td>[1, 10<sup>3</sup>]</td></tr><tr><td>l<sub>1</sub></td><td>[10<sup>-4</sup>, 1]</td></tr><tr><td>(l<sub>i</sub>)<sub>2&le;i&le;5</sub></td><td>[10<sup>-2</sup>, 10<sup>3</sup>]</td></tr><tr><td>N&#x0302;<sub>tot</sub></td><td>[10<sup>3</sup>, 10<sup>7</sup>]</td></tr><tr><td>K<sub>1</sub></td><td>[10<sup>-5</sup>, 1]</td></tr><tr><td>(&#x0303;ρ<sub>i</sub>, ρ<sub>i</sub>)<sub>1&le;i&le;4</sub></td><td>[0, 1]</td></tr><tr><td>r&#x0302;<sub>init</sub></td><td>[0, 1]</td></tr><tr><td>f&#x0302;<sub>1</sub></td><td>[0, 1]</td></tr><tr><td>n</td><td>[0, 10]</td></tr><tr><td>τ</td><td>[4, 30]</td></tr></tbody></table>

visually from Figure 3.1). Also, the half-life of the growing follicles emanating from the reserve is expected to be comparable with or little faster than that of quiescent follicles, so that we choose $l_i$, $2 \le i \le 5$ ranging from $10^{-2}$ to $10^3$.

3. The activation rate $\lambda_1$ is controlled by the most AMH-secreting follicles in compartments III and IV. We thus have to fix the (abscissa of) the inflection point of $\lambda_1$ ($1/K_1$, see (3.2)) in a range where the value of $y_3^{tot}+y_4^{tot}$ will be operative on the activation rate, namely $K_1 \sim 10^{-2}$. The Hill exponent, $n$, controls the steepness of the activation rate, and is searched in the interval $[0, 10]$. There is a 23 day-long delay needed to observe the first growing follicles emanating from the reserve after birth (Dataset C, see right panel of Figure 3.2). We thus constrain the delay, $\tau$, to be no longer than 30 days.

4. The germ cells are committed either to enter the follicle stages or to undergo massive apoptosis [58] in early life. Naked germ cells are thus depleted very quickly, so that $l_0$ (the normalized decay rate of germ cell destined to the reserve of quiescent follicles) is lower bounded by 1.

5. To fix the search interval for $\hat{N}_{tot}$, we use the biologically available information related, on one hand, to the total germ cell number at initial time, $N_{tot} = y_0(t_0) + \bar{y}_0(t_0)$, and, on the other hand, to the duration needed to build up the pool of quiescent (reserve) follicles, which can be approximated by $1/m_0$:

* At initial time $t_0$, 3 days before birth (which corresponds to embryonic day 16.5), the number of germ cells peaks around 15000 cells [27].

* The building up of the follicle reserve occurs during the first week of life in the mouse [50, 12]. This duration can be assessed accurately by following the ingression, from the ovarian epithelium, of somatic germ cells that will assemble with the germ cells to form the follicles. The ingression process is completed at 7 days of age.

Consistently, we expect $m_0$ to be on the order of $10^{-1}$ and $N_{tot}$ on the order of $10^4$. Based on these considerations, and given the uncertainty affecting the assessment of $\bar{m}_0$, we kept a rather wide search interval for $\hat{N}_{tot}$, from $10^3$ to $10^7$.
## 4.2 Error models for the three datasets

We now introduce the regression (error) model for each dataset, and define the composite cost function to be optimized.

For any dataset, the regression model can be written as :

$$S_i(t_j) = f_i(\bar{y}(t_j), y(t_j), \theta) + \mathcal{N}(0, \sigma_i(t_j)^2) \quad (3.7)$$

where,

* $t_j$ is the observation time (age),

* $S_i$ is an output variable to be compared to an experimental observation,

* $f_i$ is a transformation applied to some/all observations and the corresponding output variables,

* $\sigma_i(t)^2$ is the (empirical or estimated) variance,

* $\theta$ is the parameter vector defined by (3.6).

### 4.2.1 Dataset A: whole life data on follicle numbers

We can observe on Figure 3.1 that, in each compartment $i$ ($1 \le i \le 5$), the variance of Dataset A changes with time. There is a trend toward a larger dispersion for larger population values, which leads us to choose a statistical model with multiplicative noise. In addition, follicle numbers between compartments vary by several orders of magnitude, so that a natural transformation is logarithmic scaling. In the log-scale, the temporal trend in the dispersion seems to disappear.

To check statistically the time-independence of the variance, we apply the Breusch-Pagan test to each compartment (function het_breuschpagan from the Python package statsmodels). The null hypothesis corresponds to the homoscedastic case. First, we interpolate the discrete, log-transformed data with a spline function. Then we compute the residuals as the difference between the (log-transformed) data and the spline value at all observation times. The resulting p-values of the $\chi^2$ test ($3.95 \times 10^{-6}$; 0.28; 0.14; $8.93 \times 10^{-4}$; $2.79 \times 10^{-5}$) allow us to reject the alternative hypothesis for compartments 1 and 2, so that the homoscedasticity hypothesis is valid. For compartments 3, 4, and 5, this hypothesis appears to be not statistically valid. Since Dataset A is mostly composed of datapoints before 200 days, and the control term in the activation rate is mainly active in this time interval (thereafter the number of follicles, including the AMH secreting follicles become really low), we also performed a Breusch-Pagan restricted on the interval $[t_0, 200]$ for each compartment. All resulting p-values (0.47, 0.058, 0.083, 0.4, 0.058) are compatible with the homoscedasticity hypothesis.

We will thus stick to the following time-independent model :

$$S_i^A(t_j) = \log(y_i^{tot}(t_j)) + \mathcal{N}(0, \sigma_i^2) \quad (3.8)$$

where $\sigma_i$ is a constant parameter to be estimated.
<figure><img src="image_33.png" /><figcaption>FIGURE 3.4 – Logarithmic transformation of follicle numbers in each class as a function of age and spline-based interpolation.</figcaption></figure>

Blue dots: Data points of Dataset A. Orange solid lines: spline-based estimation of the data mean.

The green dashed lines represent the estimated standard deviation. The spline function is of degree 4 and is determined with the function *UnivariateSpline* from the Python package *scipy*.

### 4.2.2 Error model for point data comparing controlled and uncontrolled follicle activation (Dataset B)

Let $(\bar{z}_i, z_i)_{0 \le i \le 5}$ be the same model outputs as $(\bar{y}_i, y_i)_{0 \le i \le 5}$ in the case when the control term is inactive in the formulation of the activation rate ($K_1 = 0$ in (3.2)).

As a consequence, at the observation times, $t_j \in \{25, 120\}$, we deal simultaneously with data in the control and uncontrolled ($K_1 = 0$) situation:

$$ \left\{ \begin{array}{l} O_1^{tot}(t_j) = y_1^{tot}(t_j), \qquad Q_1^{tot}(t_j) = z_1^{tot}(t_j), \\ O_G(t_j) = \sum_{i=2}^{5} y_i^{tot}(t_j), \qquad Q_G(t_j) = \sum_{i=2}^{5} z_i^{tot}(t_j). \end{array} \right. $$

On the other hand, we only have the empirical mean and variance available, as the individual data points are not provided in [19]. In consistency with the regression model (3.8), we also choose a multiplicative, lognormal model, yet with known variance, that is, for $t_j \in \{25, 120\}$:

$$ \left\{ \begin{array}{l} S_1^B(t_j) = \log(y_1^{tot}(t_j)) + \mathcal{N}(0, v_1(t_j)) \\ S_2^B(t_j) = \log\left(\sum_{i=2}^{5} y_i^{tot}(t_j)\right) + \mathcal{N}(0, v_2(t_j)) \\ S_3^B(t_j) = \log(z_1^{tot}(t_j)) + \mathcal{N}(0, v_3(t_j)) \\ S_4^B(t_j) = \log\left(\sum_{i=2}^{5} z_i^{tot}(t_j)\right) + \mathcal{N}(0, v_4(t_j)) \end{array} \right. $$
where the variances $v_i(t)$ are the empirical variances.

### 4.2.3 Error model for early life data on the contribution of the first-wave follicle pool (Dataset C)

Dataset C provides the proportion of first wave follicles among the whole growing follicles, at 7 different time points:

$$P_0(t_j) = 100 \times \frac{\sum_{i=2}^{5} \bar{y}_i(t_j)}{\sum_{i=2}^{5} y_i^{tot}(t_j)}, \quad \text{with } t_j \in \{13, 23, 35, 45, 60, 90, 105\}.$$

As for dataset B, we only have the empirical mean and variance of the data. In the absence of further information on the proportion $P_0$, we choose the simplest regression model, that is an additive Gaussian model:

$$S_1^C(t_j) = P_0(t_j) + \mathcal{N}(0, v_5(t_j))$$

where the variance $v_5$ is the empirical variance except for the first and the last time points, where there is a single value ($P_0$ equals 0% or 100%), and for which we set $v_5$ to an arbitrarily low value, $v_5 = 0.01$.

### 4.2.4 Cost function

For the general regression model (3.7), with a set of observable data O, the log-likelihood function is:

$$l(S, O, \theta) := -\log \mathcal{L}(O, \theta) = \sum_{t_j} \left[ \frac{1}{2} \left( \frac{S(t_j) - O(t_j)}{\sigma_i} \right)^2 + \log(\sqrt{2\pi}\sigma_i) \right]$$

Dividing the residuals by the variance in the log-likelihood function makes it independent of the order of magnitude of the data in the O dataset. As a result, the values of the log-likelihood functions are similar between the datasets, and we combine them into the single composite cost function

$$J(\theta) = J_A(\theta) + J_B(\theta) + J_C(\theta)$$

with

$$J_A(\theta) = \sum_{i=1}^{5} l(S_i^A, O_i^{tot}, \theta)$$

$$J_B(\theta) = l(S_1^B, O_1^{tot}, \theta) + l(S_2^B, O_2, \theta) + l(S_3^B, O_3^{tot}, \theta) + l(S_4^B, O_4, \theta)$$

$$J_C(\theta) = l(S_1^C, P_0, \theta)$$
<figure><img src="image_34.png" /><figcaption>FIGURE 3.5 – Fitted follicle numbers in each compartment.</figcaption></figure>

The orange lines represent the total number of follicles in compartments I to V (from left to right and top to bottom), $y_i^{tot}$, computed with the best fit parameters. The grey areas delimit the 95 % confidence interval computed from the model error (3.8). The green lines represent the follicle numbers in the uncontrolled (AMH-KO-like) situation ($K_1 = 0$).

Rightmost bottom panel : Activation rate as function of age (red line) and number of AMH secreting follicles $y_3^{tot} + y_4^{tot}$ (blue line).

## 4.3 Parameter estimation results

### 4.3.1 Best fit parameter values

We have used the Data2Dynamics Matlab package [57] to minimize the cost function $J(\theta)$ and estimate the optimal parameter vector $\hat{\theta}$ (Table 3.3).

The best fit is illustrated on Figure 3.5 for Dataset A and Figure 3.6 for Datasets B and C. There is a general good agreement between the experimental data and estimated values, especially in a context of heterogenous data. The estimated total follicle number in each compartment is consistent with Dataset A, and most of the data-points belong to the 95% confidence interval as defined by $\sigma_i$ in (3.8). For Dataset B, the estimated values fall within or close to the 95 % confidence interval of the experimental data (red bars), except for the number of growing follicles at age 25. The estimated values capture well the most salient feature of this Dataset, i.e. the more pronounced depletion in the reserve of primordial follicles at 120 days in the AMH-KO situation compared to the wild-type. For Dataset C, the global pattern of the proportion of first wave follicles is well captured, even if the estimated curve is outside the 95 % confidence interval between 40 and 90 days.

The changes in the activation rate with age are illustrated on the rightmost bottom panel of Figure 3.5 (red curve). The first part (horizontal segment) corresponds to the delay $\tau$. At activation onset, the number of AMH secreting follicles (blue curve) is low and their feedback is not operative, so that activation occurs at maximal rate (vertical segment ending up at $7.17 \ 10^{-3} d^{-1}$). At the peak follicle number (around day 100), activation is significantly slowed down (minimum at $3.54 \ 10^{-3} d^{-1}$), and the activation rate decreases (resp. increases) from one side of the peak to the other. Finally, the follicle number falls below a threshold and the activation rate goes back to the maximum value ($7.56 \ 10^{-3} d^{-1}$). On the whole lifespan, the time average of the activation rate is $6.57 \ 10^{-3} d^{-1}$. Comparing with the AMH-KO-like situation, the feedback on the activation rate results in a gain of 200 quiescent follicles around the peak of AMH secreting follicle number.
<figure><img src="image_35.png" /><figcaption>FIGURE 3.6 – Fitting results compared with Datasets *B* and *C*.</figcaption></figure>

Left panel : The black (resp. blue) bars represent the experimental values of Dataset *B* (the same as in Figure 3.2) in the wild-type (resp. AMH-KO) situation, with the standard deviation in red. The orange and green bars correspond to the best fit values, respectively in the controlled and uncontrolled situation. Right panel : The orange line represents the estimated proportion of first-wave follicles amongst all growing follicles, superimposed on the experimental data of Dataset *C* (blue points : experimental data, red bars : standard deviation).
<table><tr><td rowspan="2">1<sup>st</sup> wave</td><td>l<sub>0</sub></td><td>l<sub>1</sub></td><td>l<sub>2</sub></td><td>l<sub>3</sub></td><td>l<sub>4</sub></td><td>l<sub>5</sub></td></tr><tr><td>3.54<br/>∈ [2.35, 5.83]</td><td>6.91<br/>∈ [5.39, 8.75]</td><td>4.79<br/>∈ [3.55, 6.38]</td><td>2.87<br/>∈ [2, 4.12]</td><td>1.61<br/>∈ [1, 4.24]</td><td>1.42<br/>∈ [1, 4.44]</td></tr><tr><td>Reserve</td><td>l<sub>0</sub></td><td>l<sub>1</sub></td><td>l<sub>2</sub></td><td>l<sub>3</sub></td><td>l<sub>4</sub></td><td>l<sub>5</sub></td></tr><tr><td></td><td>1.95<br/>∈ [1, 3.07]</td><td>0.0197<br/>∈ [0.0176, 0.0242]</td><td>1.42<br/>∈ [1.08, 1.97]</td><td>2.76<br/>∈ [2.03, 3.92]</td><td>14.3<br/>∈ [7.50, 21.5]</td><td>10.1<br/>∈ [3.25, 14.4]</td></tr></table>

<table><tr><td rowspan="2">1<sup>st</sup> wave</td><td>ρ̄<sub>1</sub></td><td>ρ̄<sub>2</sub></td><td>ρ̄<sub>3</sub></td><td>ρ̄<sub>4</sub></td></tr><tr><td>0.0359<br/>∈ [0.0234, 0.0574]</td><td>0.983<br/>∈ [0.644, 1]</td><td>0.771<br/>∈ [0.468, 1]</td><td>0.985<br/>∈ [0.388, 1]</td></tr><tr><td>Reserve</td><td>ρ<sub>1</sub></td><td>ρ<sub>2</sub></td><td>ρ<sub>3</sub></td><td>ρ<sub>4</sub></td></tr><tr><td></td><td>0.384<br/>∈ [0.267, 0.539]</td><td>1<br/>∈ [0.823, 1]</td><td>1<br/>∈ [0.490, 1]</td><td>1<br/>∈ [0.290, 1]</td></tr></table>

<table><caption>TABLE 3.3 – Best fit parameter values $\hat{\theta}$ (all parameters are expressed in linear scale) together with their confidence intervals (CI). The confidence intervals reaching the bounds are highlighted in red.</caption><tr><td>Estimated variances</td><td>σ<sub>1</sub> = 0.315<br/>∈ [0.280, 0.357]</td><td>σ<sub>2</sub> = 0.348<br/>∈ [0.310, 0.394]</td><td>σ<sub>3</sub> = 0.314<br/>∈ [0.280, 0.356]</td><td>σ<sub>4</sub> = 0.480<br/>∈ [0.425, 0.549]</td><td>σ<sub>5</sub> = 0.325<br/>∈ [0.289, 0.370]</td></tr><tr><td>Activation function</td><td>f̂<sub>1</sub> = 0.539<br/>∈ [0.307, 1]</td><td>K<sub>1</sub> = 0.00689<br/>∈ [0.00255, 0.0118]</td><td></td><td>n = 10<br/>∈ [0.713, 10]</td><td>τ = 22.3<br/>∈ [20.89, 23.99]</td></tr><tr><td>Initial conditions</td><td>N̂<sub>tot</sub> = 115 878<br/>∈ [84 333, 159 588]</td><td>τ̂<sub>init</sub> = 0.894<br/>∈ [0.808, 0.938]</td><td></td><td></td><td></td></tr></table>
4.3.2 Practical identifiability of the parameters

We have assessed the practical identifiability of the parameters from the likelihood profile ap-
proach [57]. Most of the parameters (20 parameters over 31) have a parabolic (or quasi-parabolic)
likelihood profile and are identifiable, with a clear minimum value (see the four top panels in Fi-
gure 3.7). For $(\bar{\rho}_i)_{2\le i\le 4}$, $(\rho_i)_{2\le i\le 4}$, $\bar{l}_4$, $\bar{l}_5$, $l_0$, $\tilde{f}_1$ and $n$, the likelihood profile hits one of the bounds
of the search interval with no statistically well-defined minimum (see the four bottom panels in
Figure 3.7).

We can make the following observations on the best fit parameter values together with their
confidence intervals :

* The normalized decay rates ($l_i$, $2 \le i \le 5$) of the growing follicles are clearly higher than one, which highlights a timescale difference between the dynamics of the growing follicles compared with quiescent follicles.

* The bounds of the estimated exit rate of quiescent follicles ($l_1$, 0.0176 lower bound and 0.0242 upper bound of CI) amounts to a half-life of 29-39 days ($\ln 2/l_1$), and involves a significantly higher death rate than activation rate ($\rho_1 < 0.5$ with high confidence). Given the identifiability of $l_1$, $\rho_1$, $K_1$ and the estimated lower bound on $\tilde{f}_1$, we can assess in a statistically significant way the changes in the activation rate $\lambda_1$ with aging (see bottom right panel of Figure 3.5).

Both $K_1$ and $\tilde{f}_1$ display a likelihood profile with at least two local minima, which is consistent with their sorting as locally identifiable parameters from the theoretical identifiability analysis (see Sous-section 3.2). The constraint on the upper bound of the Hill exponent, $n$, is saturated, suggesting that an ever steeper, hence heaviside-like, activation function is expected. This would suggest that there are mostly two operating modes according to the AMH-secreting follicle number : a maximal activation rate under a threshold number (which can be estimated by the abscissa of the inflexion point, $(1/\hat{K}_1 = 145)$) and a minimal rate above this threshold.

* The best fit values of $\rho_2$, $\rho_3$, $\rho_4$ are equal to 1 (and close to 1 for $\bar{\rho}_2$, $\bar{\rho}_3$, $\bar{\rho}_4$), i.e. the death rates of the growing follicles in compartments II, III and IV are negligible. This observation is consistent with the estimation of death rates performed in [23] in a piecewise linear framework, where the switching time was set to day 20, a value close to the best fit value of the delay ($\hat{\tau} = 22.3$). Setting the corresponding death rates to zero ($g_2 = g_3 = g_4 = 0$), we can compute a predicted average duration of follicle growth of about 58 days ((1/1.42+1/2.76+1/14.3)/0.0197), which is consistent with the minimal duration (40 days) observed in [66].

* Compared to the first-wave follicles, the dynamics of the growing follicles emanating from the reserve is slower, as illustrated by the best fit values of ($\bar{l}_i$, $2 \le i \le 5$) (greater than one), even if the confidence intervals of ($\bar{l}_4$) and ($\bar{l}_5$) include one. This timescale difference is consistent with the comparative assessment of the minimal transit times available in [66]. In a similar way, it is expected that $\bar{l}_1$ be far much greater than $l_1$, since the first wave follicles do not undergo a quiescent state, and that $\bar{l}_0$ be greater than $l_0$, since first wave follicles are formed approximately one week before the reserve follicles.
<figure><img src="image_36.png" /><figcaption>**FIGURE 3.7** – Representative likelihood estimation profiles parameters.
Top panels : Instances of identifiable parameters : (quasi) parabolic likelihood profiles.
Bottom panels : Instances of not practically identifiable parameters : hitting bound profiles.
The dashed red lines delimit the 95 % confidence interval.</figcaption></figure>

**4.4 Statistical model selection**

In this subsection, we perform a serie of statistical comparisons between nested models, from
simplified versions of the full model (system (3.1)). The objective is to quantify the trade-off bet-
ween the goodness of the fit and simplicity of the model (number of parameters), and to check
whether describing the dynamics of germ cells is statistically relevant.

More specifically, we define the two following sub-models. Let :

• $\mathcal{M}_1$ be the nested model with $\tau = 0$ (no delay in the activation term (see (3.2)),

• $\mathcal{M}_2$ be the nested model with $m_0 = r_0 = g_0 = 0$, $\bar{m}_0 = \bar{r}_0 = \bar{g}_0 = 0$, (no germ cell compartment, and initial conditions given by $y_1(t_0) = y_{init}$ and $\bar{y}_1(t_0) = \bar{y}_{init}$).

To compare $\mathcal{M}_1$ and $\mathcal{M}_2$ with the full model, we compute

• the Akaike information criterion (AIC) : $= 2k - 2 \log \mathcal{L}(\theta)$ with $k$ the number of estimated parameters,

• the Bayesian information criterion (BIC) : $= k \ln N - 2 \log \mathcal{L}(\theta)$ with $N = 707$ the number of observations.

The model selection is quantified by the AIC difference $\Delta_i^{\text{AIC}} = \text{AIC}_i - \min_i \text{AIC}_i$, and
the BIC difference $\Delta_i^{\text{BIC}} = \text{BIC}_i - \min_i \text{BIC}_i$, with the corresponding weights : $w_i^{\text{AIC/BIC}} = \frac{\exp(-0.5\Delta_i^{\text{AIC/BIC}})}{\sum_{i=1}^{3} \exp(-0.5\Delta_i^{\text{AIC/BIC}})}$, as well as by the p-values associated with the likelihood-ratio tests (Table 3.4).

From [9], if the “best” model corresponds to $\Delta = 0$, one can apply the following rules

• $\Delta \le 2 \Rightarrow$ the tested model is as probable as the “best” model;

• $2 < \Delta \le 7 \Rightarrow$ the tested model could be an alternative to the “best” model;

• $7 < \Delta \le 10 \Rightarrow$ the tested model is less probable that the “best” model;
<table><thead><tr><th>Model<br>(parameters number)</th><th>-2 ln L</th><th>AIC</th><th>BIC</th><th>Likelihood-ratio<br>test</th></tr></thead><tbody><tr><td>Full model (S)<br>( 31 parameters)</td><td>794.15</td><td>856.15<br>w<sup>AIC</sup> ≈ 1</td><td>997.54<br>w<sup>BIC</sup> ≈ 1</td><td></td></tr><tr><td>M<sub>1</sub><br>( 30 parameters)</td><td>1086.26</td><td>1146.26<br>Δ<sup>AIC</sup> = 290.11<br>w<sup>AIC</sup> ≈ 0</td><td>1283.09<br>Δ<sup>BIC</sup> = 285.55<br>w<sup>BIC</sup> ≈ 0</td><td>292.11<br>p<sub>1</sub> ≈ 0</td></tr><tr><td>M<sub>2</sub><br>( 29 parameters)</td><td>879.56</td><td>937.56<br>Δ<sup>AIC</sup> = 81.41<br>w<sup>AIC</sup> ≈ 0</td><td>1069.8<br>Δ<sup>BIC</sup> = 72.26<br>w<sup>BIC</sup> ≈ 0</td><td>85.40<br>p<sub>2</sub> ≈ 10<sup>-16</sup></td></tr></tbody></table>

**TABLE 3.4** – Model selection criteria used for the comparison of the full model with $M_1$ and $M_2$. Differences $\Delta_i$ and weights $w_i$ associated with the AIC and BIC criterion and p-value associated with the likelihood ratio test.

• $10 < \Delta \Rightarrow$ the tested model can be discarded.

Whatever the criterion, the model comparison selects the full model (S) as the best model. It thus appears that the delay $\tau$ is useful for capturing the follicle dynamics early after birth. Indeed, when fitting model $M_1$ to the different datasets (data not shown), the proportion of first wave follicles decrease far too rapidly and the model fails to capture the information provided by Dataset C. In addition, the follicle dynamics in compartments 2 and 3 is poorly described along the first 50 days, and there is a spurios initial peak in the curve fitting compartment 3. In the same spirit, including the dynamics of germ cells help to better explain the early follicle dynamics, which compensates for adding 2 parameters.

# 5 Biological interpretation and model predictions

## 5.1 Insight into the germ cell dynamics

The results of *sous-section 4.4* have grounded the interest to include prefollicular steps into the modeling of follicle population dynamics. A practical limit to benefit from this comprehensive modeling framework is that, for sake of parameter identifiability, we had to deal with compound (rescaled) parameters ($\hat{N}_{tot}$ and $\hat{r}_{init}$), and could not get separate information on the initial conditions, on one side, and kinetic parameters, on the other side, for the germ cell compartments.

The model and estimation results can nevertheless be used to predict the partition of the initial (total) germ cell number into medulla-located germ cells, which are committed to the first-wave, and cortex-located germ cells, which are committed to the reserve.

First, from the results of the identifiability analyses (*Sous-section 3.2* and *Sous-section 4.3*), we can predict the cumulative number of follicles entering the first wave ($\hat{m}_0\hat{y}_0/\hat{l}_0$) compatible with Dataset A, which is on the order of 15000, and the cumulative number of follicles entering the reserve ($m_0y_0/l_0$), on the order of 6300.

Then, using additional information constraining two over the four parameters $m_0$, $\bar{m}_0$, $y_0$ and $\bar{y}_0$, we can predict the value of the initial partition of germ cells, $r_{init} = \frac{\bar{y}_0(t_0)}{y_0(t_0)+\bar{y}_0(t_0)}$.

More specifically, from the biological specifications enunciated in *Point 5* from *Sous-section 4.1*, we choose, as an illustration, $m_0 = 1/10$ and $N_{tot} = 15000$, and obtain a value of $r_{init} = 0.84$. This would mean that there is an imbalance between the two follicle populations,
characterized by an excess of germ cells potentially committed to the first wave, which is related to the imbalance in the predicted germ cell inflow into each follicle subpopulation.

Note however that this prediction is very dependent on the values selected for the *a priori* information. While the assessment of the duration needed to build the quiescent pool is rather accurate and consistent between different literature sources, the assessment of the cell number at age embryonic 16.5 – based on the rare studies providing information on the maximal number of germ cells reached during oogenesis and on the developmental time when this maximum is reached [44, 52] – is much more imprecise. The peak value (at which point cell proliferation remains concomitant with cell death) is subject to a high variability resulting from between strain and inter-individual differences, and to both a low time resolution which blurs the peak location, and possible inaccuracy in the assessment of embryonic age.

Whatever it be, the predicted value of $r_{init}$ could be validated by future experiments, once they be technically feasible, given that the contrasted spatial distribution of the follicle pools in the developing ovary [12, 62] should make it possible to implement lineage tracing experiments following the germ cells until they die or get embedded into a follicle.

A more straightforward way to get this result would be to monitor the total number of germ cells more densely around birth. We now illustrate as a proof-of-concept the first step towards an experimental design procedure assisted by the model. From the best-fit parameter set listed in Table 3.3, and using the additional *a priori* information on $N_{tot}$ and $m_0$, we generate a couple of synthetic data on the germ cell numbers and merge them to the genuine experimental data. We then run a new estimation round detailed in the next subsection.

## 5.2 Towards experimental design : an instance of targeted data acquisition counseling

<figure><img src="image_37.png" /><figcaption>**FIGURE 3.8** – Synthetic data generated with the best fit parameters (Table 3.3).</figcaption></figure>

Left panel : Simulated values of $y_0^{tot}(t_j)$ noised with a log-normal distribution. Right panel : Simulated values of $p_1(t_j)$ noised with a normal distribution. The red dashed lines represent the simulated (non noised) values. The grey areas represent the 95 % confidence interval due to the noise.

The following synthetic data (see Figure 3.8) were simulated on a fixed time grid, to build the additional Dataset D :

$$ \bullet O_7(t_k) = \log y_0^{tot}(t_k) + \mathcal{N}(0, 0.1^2), (t_k = t_0 + k)_{1 \le k \le 10}, $$
* $O_8(t_l) = p_1(t_l) + \mathcal{N}(0, 0.05^2)$, with $p_1 = \frac{\bar{y}_0}{y_0+\bar{y}_0}$, $(t_l = t_0 + 3l - 2)_{1\le l\le 10}$.

For the sake of simplicity, we take integer values for the $t_k$ and $t_l$ s, while the $t_j$ values in Dataset A are real. Note that $y_0^{tot}(t_0) = N_{tot}$ and $p_1(t_0) = r_{init}$.

To perform the new estimation, we update the estimation vector as

$$ \theta = [ (\bar{l}_i)_{0\le i\le 5}, (\bar{\rho}_j)_{1\le j\le 4}, (l_i)_{0\le i\le 5}, (\rho_j)_{1\le j\le 4}, \tilde{f}_1, K_1, \tau, n, N_{tot}, r_{init}, \bar{\rho}_0, \rho_0 ] $$

in which $N_{tot}$ and $r_{init}$ replace the former parameters $\hat{N}_{tot}$ and $\hat{r}_{init}$, and $\bar{\rho}_0$, $\rho_0$ are additional dimensionless parameters :

* $\bar{\rho}_0 = \frac{\bar{m}_0}{g_0+\bar{m}_0-\bar{r}_0}$,

* $\rho_0 = \frac{m_0}{g_0+m_0-r_0}$.

$\bar{\rho}_0$ and $\rho_0$ are searched in $[0, 100]$, $r_{init}$ in $[0, 1]$, and $N_{tot}$ between $10^{3.8}$ and $10^{4.5}$ in $\log_{10}$ scale.

As illustrated on Figure 3.9, the fitting results obtained after either estimation from the initial datasets or the augmented datasets (adding Dataset D) are very similar. In addition, introducing the synthetic data enables a significant gain in parameter identifiability. In particular, all 4 parameters related to the germ cells, $N_{tot}$, $r_{init}$, $\rho_0$ and $\bar{\rho}_0$ are well estimated : the reference values fall within the estimated confidence intervals (see for example the two rightmost bottom panels of Figure 3.10). The likelihood profile of the other parameters is almost unchanged (compare panels of Figure 3.7 with those of Figure 3.10).

## 5.3 In silico experiments

To illustrate the model potential for testing functional hypotheses and perform in silico experiments, we mimic an experimental study which has investigated the effect of long term administration of AMH in mice [38, 37]. The translational application for medicine would be both related to the preservation of the follicle reserve in a physiological or iatrogenic situation (cancer chemotherapy), reversible contraception or potentiation of ovarian stimulation protocols. The experiment is based on the several-week long administration of high doses of AMH, followed by the monitoring of quiescent follicles and growing follicles up to fifty days after stopping AMH treatment. The corresponding numerical experiment is illustrated on Figure 3.11, where the period of AMH administration is materialized by the grey area. The simulation captures well the three main steps of the experiment :

1. During exposure to high levels of AMH, the activation rate is blocked, the growing compartments are not fueled anymore by the reserve and they progressively get empty. The emptying duration gives kinetic information on the duration of follicle growth [38]. In parallel, the reserve is preserved from the activation-induced depletion, so that the follicle numbers in the reserve at the end of treatment is higher in treated mice than in control mice.

2. When AMH is removed, there are no more follicles in the AMH-secreting classes, hence no feedback on the activation rate, which is no more controlled. The activation rate is maximal, as well as the rate of reserve depletion. This rapid depletion compensates for the previous preservation of the reserve, so that the follicle numbers in the reserve become comparable or even smaller than in the control case. This step is analogous to the first wave in the neonatal period, yet in an adult background.

3. As new follicles enter the AMH-secreting classes, they exert a feedback on the activation rate, which slows down, getting back to a physiological situation.
<figure><img src="image_38.png" /><figcaption>**FIGURE 3.9** – Estimation performed with the original or augmented datasets.

Top panels: Fitting follicle numbers after estimation on the initial datasets (Datasets A, B and C, orange lines) and with the augmented datasets (Datasets A, B C, and D, dotted black lines). The blue points correspond to the experimental Dataset A, the purple points to the additional Dataset D. Bottom panel: Change in ratio *p*₁ along the first 30 days; orange lines: initial datasets, black dotted lines: augmented datasets.</figcaption></figure>

The scenarios represented by the blue and green lines on Figure 3.11 differ by the preservation ratio at the end of treatment (relative number of quiescent follicles in treated versus control mice). In the first case (green scenario), the ratio is set exactly as on Figure 1B in [37], while it is a little lower in the second case (blue scenario). In the latter case, the number of quiescent follicles in the treated case (blue line) goes below the control number (orange line), while in the former the curves get closer without intersecting.

Note that both situations are realistic since (i) the experimental data are not truly kinetic; the comparison between follicle numbers at day 0 or day 30 is not performed on the same mice as follicle counting is an invasive procedure, and, in addition, only 3 replicates are available in each condition (ii) the long-term exposure to high AMH levels might have disturbed the AMH signaling
<figure><img src="image_39.png" /><figcaption>**FIGURE 3.10** – Likelihood profiles of some parameters for the estimation with synthetic data. The dashed orange lines point to the reference values. The dashed red lines delimit the 95 % confidence interval.</figcaption></figure>

system and perturbed both sensitivity to AMH (via alteration in AMH receptor expression) and
ability to secrete AMH.

**6 Discussion**

Given the major ethical and technical difficulties to number germ cells in women (except in large follicles, where the cavity can be visualized through ultrasonography [1, 3]), quantitative data on whole germ cell dynamics are extremely scarce. One therefore requires animal models to further develop fundamental knowledge in female reproductive biology, notably regarding the aging of the ovarian function. Exhaustive germ cell /ovarian follicle counting remains an invasive and extremely tedious task, which explains that densely sampled age-dependent germ cell counts are not available even in animal models, with the notable exception of the mouse. Indeed, long term and regular counting of ovarian follicles has been performed in this species [34, 23, 25]. Moreover, in the mouse, genetic manipulations have proven very useful to gain functional knowledge on the dynamics and control of follicle development. With the advent of new imaging techniques, especially those involving transparization techniques, even more data should become available, possibly in the mouse [26], but also in other vertebrate species such as fishes [41].

In this work, we have taken advantage of the existing data providing us with the follicle numbers according to different maturities on a whole life horizon. We could in particular merge different datasets harvested in the same mouse strain along different studies, which considerably enriched the data available in (very) young mice. Note that there are great differences between mouse strains [10] both in the size of the follicle reserve and in follicle dynamics, which hampers merging data from mice with different genetic backgrounds.

Despite their interest, these datasets have some limitations, which we had to cope with for our modeling and data-fitting purposes. Some are extrinsic, as the loss of the “true dynamics”, since the measurements are invasive, and measures at different ages correspond to different animals. Others are more intrinsic and correspond to the choices performed by the experimentalists and data-scientists who first handled the raw data. At a given time (hence age) the ovaries contain
<figure><img src="image_40.png" /><figcaption>**FIGURE 3.11** – In silico experiment of long-term AMH administration inspired from [38, 37]. Upper panels : changes in the follicle (or germ cell numbers) in the reference situation (orange curves) or after an episode of administration of high dose AMH between days 80 and 120 (grey rectangles). Green lines : scenario with a high preservation rate of the reserve (1855 at day 120 to be compared to 1408); blue lines : scenario with a moderate preservation rate of the reserve (1500 at day 120). Lower panel : zoom on the dynamics of quiescent follicles. Note that the blue curve intersects the orange curve around day 140.</figcaption></figure>

follicles in different maturation stages. This implies that some of the data are paired; different datapoints should be associated with a same individual. Unfortunately, information on the pairing is lost in the current data. In addition, the authors have customized the usual follicle classification to sort their data. From the original classes V and VI of [45], they introduced a new class V pooling antral follicles with non antral follicles. The authors' motivation was to avoid possible variations induced by the ovarian cycles. There are indeed cyclic changes in the number of large antral follicles [55], yet they do not affect the smallest ones, which are far more numerous. In contrast, in mouse the dynamics of antral follicles is very contrasted with respect to the other
stages; they are in the so-called terminal phase of follicle development, whose control differs radically from the previous, so-called basal phase of development. Indeed, antral follicles becomes totally dependent on the supply of pituitary reproductive hormones, and they interact with each other through complex endocrine feedback loops within the gonadotropic axis. Not being able of handling separately antral follicles is penalizing, as most functional studies, including those used for our additional specifications [19, 66], separate them from the other follicles.

Given the nature of the data, our model focuses on activation and basal development, and we have paid a special attention to account for up-to-date knowledge of activation control. A proper use of the follicle reserve is indeed primordial for preserving fertility a sufficient long time, and it is now well-known that AMH secreted from growing follicles has a prominent role in avoiding precocious depletion of the reserve. This has not only been established through studies on the AMH-KO in mice, but also on cultures of ovarian cortex supplemented or not with AMH [18], and in cutting edge studies with transient overexpression of AMH [38]. AMH effects are not restricted to the activation step. AMH also modulates the growth speed in growing follicles and moderates the differentiation rate in mature follicles. Due to the scarcity of data, we did not introduce controlled growth rates between the growing compartments. Another way of refining the model would be to weight the contributions of the different classes in the feedback terms entering the activation rate, but, again, the sparse data available are not compatible with adding other parameters. Note also that if the number of follicle cells is higher in class IV follicles than in class III, the per cell contribution is lower [51, 15], which roughly compensates for the resulting AMH output, and substantiates the absence of weights in our formulation.

We have extended the initial modeling framework introduced by JM Faddy and co-authors to better characterize follicle development very early in life, from the postnatal period to several weeks after birth. The extended model takes into account the existence of two different populations of follicles at birth, and also includes an additional pre-follicular compartment, which allows us to account in a compact way for events occurring before birth [27]. J.M. Faddy and co-authors had to cope in a phenomenological way with the strong differences affecting follicle development in very young mice compared to elder ones, since they chose to fit piecewise linear models and introduced a switching time (at day 20 for the datasets selected for our study) to fit separately the initial part of the data. The choice of the switching time was empirical; it was chosen as the nadir of a spline-based regression curve applied to class II follicles [24].

Altogether, our knowledge-augmented approach for data-fitting fits successfully the experimental data from all three input datasets, while accounting for the most recent biological knowledge. Making use of the AMH-KO data [19] guarantees us to select parameter regimes in which the control of the activation rate is active (see the graph of $\lambda_1(t)$ as a function of time or class III and IV follicle numbers, see Figure 3.5). In addition, accounting for a controlled activation rate allows us to mimic the effects of AMH deprivation, in a permanent way (such as the AMH-KO, see the green curves on Figure 3.5) or transitory way (induction of a first-wave like dynamics in adulthood after a prolonged exposure to high AMH levels, see Figure 3.11).

We have also revealed that accounting for the dynamics of germ cells in pre-follicular stages help to monitor finely the follicle dynamics, especially in the earliest stages, occurring around birth in the mouse (and during embryonic life in larger mammal species including humans), and to describe separately follicles emanating either from the first wave or the reserve. Even if the objective of the current work was not to model the formation of follicles, we have shown that model-assisted acquisition of well targeted data would be sufficient to quantify the transition from the pre-follicular to the follicular stages. The oogenesis process leading to a follicle-embedded
oocyte from the very primordial germ cells would also benefit from a population dynamics approach
and will be the matter of further modeling work.

# 7 Acknowledgements

We wish to thank Gleb Pogudin for helpful discussions and counsels on the theoretical identi-
fiability analysis. FC and RY acknowledge support from ANSES, PNR EST 2020, project GINFIZ
(Anses 2020/01/133), as well as from INRAE, project IMMO (metaprogram DIGIT-BIO) and ANR
project OVOPAUSE (AAPG 2022).
<figure><img src="image_41.png" /></figure>
# Discussion et perspectives

Dans cette thèse, nous avons introduit un modèle de population structurée par l’âge afin de représenter le développement folliculaire tout au long de la vie reproductive. En tenant compte des différences d’abondance et de vitesse de maturation au sein de cette population, nous avons séparé ce modèle selon deux échelles de temps ce qui nous a permis de prouver l’existence et l’unicité d’un modèle limite (Chapitre 1). Ce modèle limite consiste en une EDO caractérisant l’épuisement progressif de la réserve de follicules primordiaux au cours du temps, couplée avec la distribution stationnaire d’une chaîne de Markov représentant la répartition des follicules dans les différents stades de croissance. Le couplage entre ces deux objets a lieu à travers la moyennisation de la distribution (la partie rapide du modèle) qui intervient dans l’équation différentielle régissant la dynamique de la réserve (la partie lente).

D'un point de vue numérique, la suite du chapitre 1 détaille deux algorithmes qui ont été envisagés pour simuler le modèle limite obtenu avec deux échelles de temps. À l'aide de l'algorithme le plus adapté à notre modèle, nous avons illustré numériquement la convergence ainsi que certaines grandeurs obtenues grâce au modèle qui peuvent être interprétées biologiquement.

Nous avons ensuite complexifié le modèle initial en ajoutant une troisième échelle de temps (Chapitre 2). Cet ajout permet de mieux prendre en compte l'ensemble de la dynamique folliculaire : l'échelle de temps la plus lente rend compte de l'épuisement de la réserve, l'échelle intermédiaire de la maturation des follicules au cours de la croissance basale, et l'échelle rapide de l'accélération de la maturation des follicules lorsqu'ils atteignent les stades terminaux. En utilisant les mêmes outils mathématiques que dans le chapitre 1, nous avons pu montrer l'existence et l'unicité d'un modèle limite tenant compte de la dynamique de chaque échelle de temps. Ce modèle limite est constitué de trois parties couplées entre elles (une pour chaque échelle de temps) : une EDO 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 le même esprit que pour le modèle à deux échelles de temps, les deux échelles de temps les plus rapides interviennent sur l'échelle la plus lente via une moyennisation.

Pour réussir à démontrer ce résultat, nous avons utilisé des hypothèses plus fortes sur la monotonie des taux de transition du modèle. Celles-ci assurent l'existence et l'unicité du point fixe de lu système d'EDOs caractérisant l'échelle intermédiaire. Dans un cadre général (dimension et taux de transition quelconques), l'étude des points fixes d'une EDO est particulièrement ardue ce qui explique ces hypothèses. Par contre, lorsque la dimension et la forme des taux sont fixés, il est possible qu'il ne soit pas nécessaire d'avoir des hypothèses si contraignantes. De plus, l'existence d'un point fixe peut être vérifiée au moins numériquement.

Avant de pouvoir estimer les paramètres de notre modèle, nous avons du chercher des jeux de
données compatibles entre eux malgré leurs natures hétérogènes. Nous avons alors adapté notre modèle de population structurée aux 4 jeux de données sur la souris retenus (Chapitre 3). Le modèle ainsi obtenu tient compte de l'existence d'une première vague de follicules (permettant de mettre en place la régulation) et capture plus finement le début du processus de formation des follicules (en ajoutant un compartiment en amont de la dynamique folliculaire). De plus, devant le peu de données disponibles, nous avons préféré réaliser l'estimation sur une version déterministe du modèle (système d'EDOs). En tenant compte des résultats d'une étude de l'identifiabilité théorique des paramètres, nous avons reparamétré le modèle et choisi une stratégie d'estimation (choix des modèles d'erreur et de la fonction de coût). Nous avons alors étudié l'identifiabilité théorique et pratique des paramètres du modèle.

D'un point de vue biologique, cette étude nous a permis de prédire une grandeur biologique en tenant compte de connaissances quantitatives sur la folliculogénèse chez la souris. Enfin, le potentiel du modèle a été illustré par la reproduction d'expériences biologiques *in silico*.

Dans cette thèse, nous avons considéré des modèles de population structurée par une variable discrète correspondant à des classifications biologiques des follicules. Il existe différentes classifications des follicules basées sur des critères morphologiques (par exemple [56, 30]). Afin de s'affranchir des classes du développement folliculaire, le modèle pourrait être raffiné en rendant la variable structurante continue. L'existence d'un modèle limite serait alors une question intéressante faisant intervenir des outils différents. D'un point de vue biologique, cette variable structurante continue pourrait être le diamètre des follicules. Dans ce nouveau cadre, il devient par contre difficile de définir correctement les interactions entre les follicules car celles-ci dépendent plus de critères morphologiques que du diamètre des follicules.

Du point de vue de la reproduction, l'ensemble du développement folliculaire est important mais les stades terminaux et plus particulièrement le cycle ovarien ont une dynamique sensiblement différente faisant intervenir des objets biologiques de nature différente (corps jaune, hormones...). Une première étape pour ajouter la dynamique du cycle ovarien serait d'adapter le modèle à trois échelles de temps (chapitre 2) afin que l'échelle de temps la plus rapide représente au mieux le cycle ovarien. La principale difficulté à surmonter est le choix de fonctions de taux permettant de reproduire le caractère cyclique. D'un point de vue biologique, de nombreuses et diverses études ont été réalisées sur le cycle ovarien [43, 47, 55] ; permettant d'envisager dans un second temps une estimation des paramètres de ce nouveau modèle.
Bibliographie

[1] G. Adams, J. Singh, and A. Baerwald. Large animal models for the study of ovarian follicular dynamics in women. *Theriogenology*, 78(8):1733-1748, Nov. 2012.

[2] J. Azema, M. Kaplan-Duflo, and D. Revuz. Mesure invariante sur les classes récurrentes des processus de Markov. *Z. Wahrscheinlichkeitstheorie verw Gebiete*, 8(3):157-181, 1967.

[3] A. R. Baerwald, G. P. Adams, and R. A. Pierson. Ovarian antral folliculogenesis during the human menstrual cycle : A review. *Human Reproduction Update*, 18(1):73-91, Jan. 2012.

[4] G. Ballif, F. Clément, and R. Yvinec. Averaging of a stochastic slow-fast model for population dynamics : Application to the development of ovarian follicles. *SIAM J. Appl. Math.*, 82(1):359-380, Feb. 2022.

[5] H. Baumann and W. Sandmann. Numerical solution of level dependent quasi-birth-and-death processes. *Procedia Computer Science*, 1(1):1561-1569, May 2010.

[6] C. Bonnet, K. Chahour, F. Clément, M. Postel, and R. Yvinec. Multiscale population dynamics in reproductive biology : Singular perturbation reduction in deterministic and stochastic models. *ESAIM Proc. Surveys*, 67:28, Apr. 2020.

[7] L. Bright and P. Taylor. Calculating the equilibrium distribution in level dependent quasi-birth-and-death processes. *Communications in Statistics. Stochastic Models*, 11(3):497-525, Jan. 1995.

[8] F. J. Broekmans, M. R. Soules, and B. C. Fauser. Ovarian Aging : Mechanisms and Clinical Consequences. *Endocrine Reviews*, 30(5):465-493, Aug. 2009.

[9] K. P. Burnham, D. R. Anderson, and K. P. Burnham. *Model Selection and Multimodel Inference : A Practical Information-Theoretic Approach*. Springer, New York, 2nd ed edition, 2002.

[10] J. Canning, Y. Takai, and J. L. Tilly. Evidence for Genetic Modifiers of Ovarian Follicular Endowment and Development from Studies of Five Inbred Mouse Strains. *Endocrinology*, 144(1):9-12, Jan. 2003.

[11] F. Clément and D. Monniaux. Mathematical modeling of ovarian follicle development : A population dynamics viewpoint. *Current Opinion in Endocrine and Metabolic Research*, 18:54-61, June 2021.

[12] M. H. Cordeiro, S.-Y. Kim, K. Ebbert, F. E. Duncan, J. Ramalho-Santos, and T. K. Woodruff. Geography of Follicle Formation in the Embryonic Mouse Ovary Impacts Activation Pattern During the First Wave of Folliculogenesis1. *Biology of Reproduction*, 93(4), Oct. 2015.

[13] J. E. Coxworth and K. Hawkes. Ovarian follicle loss in humans and mice : Lessons from statistical model comparison. *Hum. Reprod.*, 25(7):1796-1805, July 2010.

[14] C. Dellacherie and P. A. Meyer. *Probabilities and Potential*. Number 29, <72, 151> in North-Holland Mathematics Studies. North-Holland Publication, Elsevier, Paris, 1978.
[15] M. M. Devillers, F. Petit, V. Cluzet, C. M. François, F. Giton, G. Garrel, J. Cohen-Tannoudji, and C. J. Guigon. FSH inhibits AMH to support ovarian estradiol synthesis in infantile mice. *Journal of Endocrinology*, 240(2):215-228, Feb. 2019.

[16] R. Dong, C. Goodbrake, H. Harrington, and P. G. Differential elimination for dynamical models via projections with applications to structural identifiability. *preprint*, 2021.

[17] J. Donnez and M.-M. Dolmans. Transplantation of ovarian tissue. *Best Practice & Research Clinical Obstetrics & Gynaecology*, 28(8):1188-1197, Nov. 2014.

[18] A. L. L. Durlinger, M. J. G. Gruijters, P. Kramer, B. Karels, H. A. Ingraham, M. W. Nachtigal, J. T. J. Uilenbroek, J. A. Grootegoed, and A. P. N. Themmen. Anti-Mullerian Hormone Inhibits Initiation of Primordial Follicle Growth in the Mouse Ovary. *Hum. Reprod.*, 143(3):1076-1084, 2002.

[19] A. L. L. Durlinger, P. Kramer, B. Karels, F. H. de Jong, J. T. J. Uilenbroek, J. A. Grootegoed, and A. P. N. Themmen. Control of Primordial Follicle Recruitment by Anti-Müllerian Hormone in the Mouse Ovary1. *Endocrinology*, 140(12):5789-5796, Dec. 1999.

[20] P. Dürrenberger, A. Gupta, and M. Khammash. A finite state projection method for steady-state sensitivity analysis of stochastic reaction networks. *J. Chem. Phys.*, 150(13):134101, Apr. 2019.

[21] S. N. Ethier and T. G. Kurtz. *Markov Processes: Characterization and Convergence*. Wiley Series in Probability and Statistics. Wiley Interscience, Hoboken, 2005.

[22] M. J. Faddy and R. G. Gosden. Physiology: A mathematical model of follicle dynamics in the human ovary. *Human Reproduction*, 10(4):770-775, Apr. 1995.

[23] M. J. Faddy, R. G. Gosden, and R. G. Edwards. Ovarian follicle dynamics in mice: A comparative study of three inbred strains and an F1 hybrid. *J. Endocrinol.*, 96(1):23-33, Jan. 1983.

[24] M. J. Faddy, E. C. Jones, and R. G. Edwards. An analytical model for ovarian follicle dynamics. *J. Exp. Zool.*, 197(2):173-185, Aug. 1976.

[25] M. J. Faddy, E. Telfer, and R. G. Gosden. The kinetics of pre-antral follicle development in ovaries of CBA/Ca mice during the first 14 weeks of life. *Cell Prolif*, 20(6):551-560, Nov. 1987.

[26] Y. Feng, P. Cui, X. Lu, B. Hsueh, F. Möller Billig, L. Zarnescu Yanez, R. Tomer, D. Boerboom, P. Carmeliet, K. Deisseroth, and A. J. W. Hsueh. CLARITY reveals dynamics of ovarian follicular architecture and vasculature in three-dimensions. *Sci Rep*, 7(1):44810, Apr. 2017.

[27] J. K. Findlay, K. J. Hutt, M. Hickey, and R. A. Anderson. How Is the Number of Primordial Follicles in the Ovarian Reserve Established? *Biology of Reproduction*, 93(5), Nov. 2015.

[28] D. T. Gillespie. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. *Journal of Computational Physics*, 22(4):403-434, Dec. 1976.

[29] P. W. Glynn and S. P. Meyn. A Liapounov bound for solutions of the Poisson equation. *Ann. Probab.*, 24(2):916-931, Apr. 1996.

[30] A. Gougeon. Regulation of Ovarian Follicular Development in Primates: Facts and Hypotheses. *Endocrine Reviews*, 17(2):121-155, Apr. 1996.

[31] A. Gougeon and G. B. N. Chainy. Morphometric studies of small follicles in ovaries of women at different ages. *Reproduction*, 81(2):433-442, Nov. 1987.
[32] A. N. Hirshfield. Heterogeneity of Cell Populations that Contribute to the Formation of Pri-mordial Follicles in Rats1. *Biology of Reproduction*, 47(3):466-472, Sept. 1992.

[33] H. Hong, A. Ovchinnikov, G. Pogudin, and C. Yap. Global Identifiability of Differential Models. *Comm. Pure Appl. Math.*, 73(9):1831-1879, Sept. 2020.

[34] E. C. Jones and P. L. Krohn. The relationships between age, numbers of oocytes and fertility in virgin and multiparous mice. *J. Endocrinol.*, page 30, Feb. 1961.

[35] H.-W. Kang and T. G. Kurtz. Separation of time-scales and model reduction for stochastic reaction networks. *Ann. Appl. Probab.*, 23(2):529-583, Apr. 2013.

[36] H.-W. Kang, T. G. Kurtz, and L. Popovic. Central limit theorems and diffusion approximations for multiscale Markov chain models. *Ann. Appl. Probab.*, 24(2):721-759, Apr. 2014.

[37] M. Kano, J. Y. Hsu, H. D. Saatcioglu, N. Nagykery, L. Zhang, M. E. Morris Sabatini, P. K. Donahoe, and D. Pépin. Neoadjuvant Treatment With Müllerian-Inhibiting Substance Syn-chronizes Follicles and Enhances Superovulation Yield. *Journal of the Endocrine Society*, 3(11):2123-2134, Nov. 2019.

[38] M. Kano, A. E. Sosulski, L. Zhang, H. D. Saatcioglu, D. Wang, N. Nagykery, M. E. Sabatini, G. Gao, P. K. Donahoe, and D. Pépin. AMH/MIS as a contraceptive that protects the ovarian reserve during chemotherapy. *Proc Natl Acad Sci USA*, 114(9):E1688-E1697, Feb. 2017.

[39] J. Kuntz, P. Thomas, G.-B. Stan, and M. Barahona. Stationary distributions of continuous-time Markov chains : A review of theory and truncation-based approximations. *arXiv:1909.05794 [cond-mat, q-bio]*, Sept. 2019.

[40] T. G. Kurtz. Averaging for martingale problems and stochastic approximation. In I. Karatzas and D. Ocone, editors, *Applied Stochastic Analysis*, volume 177, pages 186-209. Springer-Verlag, Berlin/Heidelberg, 1992.

[41] M. Lesage, M. Thomas, J. Bugeon, A. Branthonne, S. Gay, E. Cardona, M. Haghebaert, F. Mahé, J. Bobe, and V. Thermes. C-ECi : A CUBIC-ECi combined clearing method for three-dimensional follicular content analysis in the fish ovary†. *Biology of Reproduction*, 103(5):1099-1109, Oct. 2020.

[42] T. Lindvall. *Lectures on the Coupling Method*. Dover Publications, Inc, Mineola, N.Y, 2002.

[43] C. Lliberos, S. H. Liew, P. Zareie, N. L. La Gruta, A. Mansell, and K. Hutt. Evaluation of inflammation and follicle depletion during ovarian ageing in mice. *Sci Rep*, 11(1):278, Dec. 2021.

[44] S. Malki, M. E. Tharp, and A. Bortvin. A Whole-Mount Approach for Accurate Quantitative and Spatial Assessment of Fetal Oocyte Dynamics in Mice1. *Biology of Reproduction*, 93(5), Nov. 2015.

[45] A. M. Mandl and S. Zuckerman. Numbers of normal and atretic oocytes in unilaterally spayed rats. *J Endocrinol*, 7(2):112-119, June 1951.

[46] S. Meleard and V. Bansaye. *Some Stochastic Models for Structured Populations : Scaling Limits and Long Time Behavior*. Springer International Publishing, 2015.

[47] H. Meunier, S. B. Cajander, V. J. Roberts, C. Rivier, P. E. Sawchenko, A. J. W. Hsueh, and W. Vale. Rapid Changes in the Expression of Inhibin α-, βA-, and βB-Subunits in Ovarian Cell Types During the Rat Estrous Cycle. *Molecular Endocrinology*, 2(12):1352-1363, Dec. 1988.
[48] S. P. Meyn and R. L. Tweedie. Stability of Markovian processes III : Foster-Lyapunov criteria for continuous-time processes. *Advances in Applied Probability*, 25(3):518-548, Sept. 1993.

[49] D. Monniaux, B. Mandon-Pépin, and P. Monget. L'atrésie folliculaire, un gaspillage programmé. *Med Sci (Paris)*, 15(2):157, 1999.

[50] L. Mork, D. M. Maatouk, J. A. McMahon, J. J. Guo, P. Zhang, A. P. McMahon, and B. Capel. Temporal Differences in Granulosa Cell Specification in the Ovary Reflect Distinct Follicle Fates in Mice1. *Biology of Reproduction*, 86(2), Feb. 2012.

[51] A. Munsterberg and R. Lovell-Badge. Expression of the mouse anti-mullerian hormone gene suggests a role in both male and female sexual differentiation. *Development*, 113(2):613-624, Oct. 1991.

[52] M. Myers, F. H. Morgan, S. H. Liew, N. Zerafa, T. U. Gamage, M. Sarraj, M. Cook, I. Kapic, A. Sutherland, C. L. Scott, A. Strasser, J. K. Findlay, J. B. Kerr, and K. J. Hutt. PUMA regulates germ cell loss and primordial follicle endowment in mice. *REPRODUCTION*, 148(2):211-219, Aug. 2014.

[53] M. F. Neuts and F. P. Kelly. Reversibility and Stochastic Networks. *Journal of the American Statistical Association*, 76(374):492, June 1981.

[54] W. Niu and A. C. Spradling. Two distinct pathways of pregranulosa cell differentiation support follicle formation in the mouse ovary. *Proc. Natl. Acad. Sci. U.S.A.*, 117(33):20015-20026, Aug. 2020.

[55] A. Numazawa and S. Kawashima. Morphometric studies on ovarian follicles and corpora lutea during the oestrous cycle in the mouse. *Reproduction*, 64(2):275-283, Mar. 1982.

[56] T. Pedersen and H. Peters. Proposal for a classification of oocytes and follicles in the mouse ovary. *Reproduction*, 17(3):555-557, Dec. 1968.

[57] A. Raue, B. Steiert, M. Schelker, C. Kreutz, T. Maiwald, H. Hass, J. Vanlier, C. Tönsing, L. Adlung, R. Engesser, W. Mader, T. Heinemann, J. Hasenauer, M. Schilling, T. Höfer, E. Klipp, F. Theis, U. Klingmüller, B. Schöberl, and J. Timmer. Data2Dynamics : A modeling environment tailored to parameter estimation in dynamical systems. *Bioinformatics*, 31(21):3558-3560, Nov. 2015.

[58] P. Rodrigues, D. Limback, L. K. McGinnis, C. E. Plancha, and D. F. Albertini. Multiple mechanisms of germ cell loss in the perinatal mouse ovary. *REPRODUCTION*, 137(4):709-720, Apr. 2009.

[59] W. Rudin. *Real and Complex Analysis*. McGraw-Hill, New York, 3rd ed edition, 1987.

[60] G. B. Schaalje and H. R. van der Vaart. Relationships Among Recent Models for Insect Population Dynamics with Variable Rates of Development. In J. Berger, S. Fienberg, J. Gani, K. Krickeberg, B. Singer, L. L. McDonald, B. F. J. Manly, J. A. Lockwood, and J. A. Logan, editors, *Estimation and Analysis of Insect Populations*, volume 55, pages 299-312. Springer New York, New York, NY, 1989.

[61] M. Strugarek, L. Dufour, N. Vauchelet, L. Almeida, B. Perthame, and D. A. M. Villela. Oscillatory regimes in a mosquito population model with larval feedback on egg hatching. *Journal of Biological Dynamics*, 13(1):269-300, Jan. 2019.

[62] H. Suzuki, M. Kanai-Azuma, and Y. Kanai. From Sex Determination to Initial Folliculogenesis in Mammalian Ovaries: Morphogenetic Waves along the Anteroposterior and Dorsoventral Axes. *Sex Dev*, 9(4):190-204, 2015.
[63] C. Thibault, M.-C. Levasseur, and R. H. F. Hunter, editors. *Reproduction in Mammals and Man*. INRA Editions, Paris, ellipses, engl. ed edition, 1993.

[64] W. H. B. Wallace and T. W. Kelsey. Human Ovarian Reserve from Conception to the Menopause. *PLoS ONE*, 5(1):e8772, Jan. 2010.

[65] E. Weinan, D. Liu, and E. Vanden-Eijnden. Nested stochastic simulation algorithms for chemical kinetic systems with multiple time scales. *Journal of Computational Physics*, 221(1):158-180, Jan. 2007.

[66] W. Zheng, H. Zhang, N. Gorre, S. Risal, Y. Shen, and K. Liu. Two classes of ovarian primordial follicles exhibit distinct developmental dynamics and physiological functions. *Human Molecular Genetics*, 23(4):920-928, Feb. 2014.

[67] X. Zhu, D. B. Finlay, M. Glass, and S. B. Duffull. Model-free and kinetic modelling approaches for characterising non-equilibrium pharmacological pathway activity : Internalisation of cannabinoid CB1 receptors. *British Journal of Pharmacology*, 176(14):2593-2607, 2019.
**Titre** : Modélisation et analyse stochastique multi-échelle en dynamique de population : application à la population folliculaire ovarienne.

**Mots clés** : modèle à compartiments, multi-échelle, chaîne de Markov à temps continu, perturbations singularières, estimation de paramètres, biologie de la reproduction

**Résumé** : Cette thèse vise à élaborer et analyser mathématiquement des modèles de populations structurées pour étudier l'évolution de la population de follicules ovariens tout au long de la vie reproducitive des mammifères. Dans un premier temps, nous avons construit un modèle de naissance, mort et migration; structuré en compartiments : chaque stade de maturation des follicules est représenté par un compartiment. Afin de tenir compte de possibles interactions, les taux de transition de cette chaîne de Markov à temps continu dépendent de l'état de la population. En utilisant une différence d'abondance et de vitesse de transit au sein des compartiments, nous obtenons un modèle dit lent/rapide. Nous démontrons alors l'existence et l'unicité d'un modèle limite dans le cadre de la théorie des perturbations stochastiques singulières en utilisant un critère de Foster-Lyapunov et une équation de Poisson bien choisie. Dans un second temps, nous ajoutons au modèle une échelle de temps intermédiaire pour modéliser plus finement l'accélération continue de la croissance des follicules. À l'aide des mêmes outils mathématiques, nous prou-

vons l'existence et l'unicité d'un modèle limite plus complexe que le précédent. Chaque échelle dépend de la précédente via une moyennisation ce qui rend l'étude générale de ce modèle particulièrement difficile. Dans le cas de 2 ou de 3 échelles de temps, nous avons ensuite poussé l'étude de la limite à l'ordre supérieur en caractérisant l'erreur entre le modèle et sa limite. Enfin, nous avons introduit un modèle (un système d'équations différentielles ordinaires) adapté aux données biologiques expérimentales disponibles chez la souris. Ce modèle capture plus finement les premières étapes de formation des follicules (par l'ajout d'un compartiment en amont de la dynamique folliculaire) et tient compte de l'existence de deux sous-populations distinctes. En tenant compte d'une étude de l'identifiabilité théorique des paramètres, nous avons reparamétré le modèle et estimé les nouveaux paramètres du modèle. Nous avons finalement illustré le potentiel du modèle de deux façons : en suggérant le type de données biologiques à ajouter afin d'identifier un plus grand nombre de paramètres et en reproduisant des expériences in silico.

**Title** : Stochastic multiscale modeling and analysis of structured population dynamics : application to the ovarian follicle population.

**Keywords** : compartmental model, multiscale, continuous time Markov chain, singular perturbations, parameter estimation, reproductive biology

**Abstract** : This thesis aims to develop and analyze mathematically structured population models to explore the dynamics of the ovarian follicle population throughout the reproductive life of mammals. Firstly, we have built a birth, migration and death stochastic process structured in compartments: each follicle maturation stage is characterized by a compartment. In order to take into account possible interactions, some transition rates of the resulting continuous-time Markov chain depend on the whole or part of the population. Using a difference in abundances and time scales within the population, we have obtained a so-called slow/fast model. We have then proved the existence and uniqueness of a limit model in the framework of stochastic singular perturbations by using a Foster-Lyapunov criterion and a well-chosen Poisson equation. In a second step, we have added to the model an intermediate timescale to refine the modeling of the continued acceleration of follicle growth. With the same mathematical tools, we have proven the existence and uniqueness of a more

complex limit model than the previous one. Each timescale depends on the previous one through an averaging, which makes the study of this model particularly difficult in a general framework. In the case of two and three timescales, we have then extended the study of the limit to the higher order by identifying the mean and variance of the error between the model and its limit. Finally, we have introduced a model (a system of ordinary-differential equation) adapted to the biological experimental data available in the mouse. This model captures more accurately the early step of follicle formation (by adding a compartment upstream the follicle dynamics) and takes into account the existence of two distinct subpopulations of follicles. Taking into consideration a study of theoretical identifiability of the parameters, we have reparameterized the model and estimated the new parameters. We have illustrated the model potential in two ways: by suggesting the type of biological data needed to identify a larger number of parameters and by reproducing in silico experiments.