Stochastic nonlinear model for somatic cell population dynamics during ovarian follicle activation

TL;DR

This paper develops a stochastic nonlinear model to analyze the dynamics of somatic cells during ovarian follicle activation, revealing the importance of cell awakening over proliferation and the impact of feedback mechanisms.

Contribution

It introduces a novel stochastic model for follicle activation, incorporating feedback effects and providing insights into the timing of cell awakening versus proliferation.

Findings

Feedback reduces proliferative cell numbers at extinction.

Model calibration using likelihood approaches.

Awakening precedes proliferation during activation.

Abstract

In mammals, female germ cells are sheltered within somatic structures called ovarian follicles, which remain in a quiescent state until they get activated, all along reproductive life. We investigate the sequence of somatic cell events occurring just after follicle activation, starting by the awakening of precursor somatic cells, and their transformation into proliferative cells. We introduce a nonlinear stochastic model accounting for the joint dynamics of the two cell types, and allowing us to investigate the potential impact of a feedback from proliferative cells onto precursor cells. To tackle the key issue of whether cell proliferation is concomitant or posterior to cell awakening, we assess both the time needed for all precursor cells to awake, and the corresponding increase in the total cell number with respect to the initial cell number. Using the probabilistic theory of first passage times, we design a numerical scheme based on a rigorous Finite State Projection and coupling techniques to compute the mean extinction time and the cell number at extinction time. We find that the feedback term clearly lowers the number of proliferative cells at the extinction time. We calibrate the model parameters using an exact likelihood approach. We carry out a comprehensive comparison between the initial model and a series of submodels, which helps to select the critical cell events taking place during activation, and suggests that awakening is prominent over proliferation.

Tables

Table 1. Table 1: Model comparison analysis. For each experimental subset and each submodel, we compute both the Akaike information criterion (AIC) and Bayesian information criterion (BIC), the AIC and BIC differences Δ i A ​ I ​ C := A ​ I ​ C i − A ​ I ​ C min assign subscript superscript Δ 𝐴 𝐼 𝐶 𝑖 𝐴 𝐼 subscript 𝐶 𝑖 𝐴 𝐼 subscript 𝐶 \Delta^{AIC}_{i}:=AIC_{i}-AIC_{\min} and Δ i B ​ I ​ C = B ​ I ​ C i − B ​ I ​ C min subscript superscript Δ 𝐵 𝐼 𝐶 𝑖 𝐵 𝐼 subscript 𝐶 𝑖 𝐵 𝐼 subscript 𝐶 \Delta^{BIC}_{i}=BIC_{i}-BIC_{\min} , and the corresponding Akaike and Bayesian weights w i A ​ I ​ C = exp ⁡ ( − 0.5 ​ Δ i A ​ I ​ C ) ∑ k = 1 6 exp ⁡ ( − 0.5 ​ Δ k A ​ I ​ C ) subscript superscript 𝑤 𝐴 𝐼 𝐶 𝑖 0.5 subscript superscript Δ 𝐴 𝐼 𝐶 𝑖 superscript subscript 𝑘 1 6 0.5 subscript superscript Δ 𝐴 𝐼 𝐶 𝑘 w^{AIC}_{i}=\frac{\exp(-0.5\Delta^{AIC}_{i})}{\sum_{k=1}^{6}\exp(-0.5\Delta^{AIC}_{k})} and w i B ​ I ​ C = exp ⁡ ( − 0.5 ​ Δ i B ​ I ​ C ) ∑ k = 1 6 exp ⁡ ( − 0.5 ​ Δ k B ​ I ​ C ) subscript superscript 𝑤 𝐵 𝐼 𝐶 𝑖 0.5 subscript superscript Δ 𝐵 𝐼 𝐶 𝑖 superscript subscript 𝑘 1 6 0.5 subscript superscript Δ 𝐵 𝐼 𝐶 𝑘 w^{BIC}_{i}=\frac{\exp(-0.5\Delta^{BIC}_{i})}{\sum_{k=1}^{6}\exp(-0.5\Delta^{BIC}_{k})} following ( burnham_model_2003, , (Chapter 2 and 3)) . The best models are highlighted in red and the remaining selected models in blue (details are provided in Appendix 6.6 ).

	Wild-Type			Mutant
Model	-$\log ℒ(\theta ;𝐱)$	AIC	BIC	-$\log ℒ(\theta ;𝐱)$	AIC	BIC
$({ℛ}_1,{ℛ}_4)$	172.87	349.74 $w=0.02$ $\Delta =7.6$	354.74 $w=0.15$ $\Delta =3.0$	149.97	303.94 $w=0.08$ $\Delta =8.8$	308.73 $w=0.03$ $\Delta =6.4$
$({ℛ}_1,{ℛ}_3)$	245.54	495.09 $\Delta >>10$	500.09 $\Delta >>10$	230.17	464.34 $\Delta >>10$	469.13 $\Delta >>10$
$({ℛ}_1,{ℛ}_2,{ℛ}_4)$	172.77	351.54 $w=0.008$ $\Delta =9.44$	359.04 $\Delta =7.4$	148.14	302.27 $w=0.02$ $\Delta =7.1$	309.46 $w=0.02$ $\Delta =7.1$
$({ℛ}_1,{ℛ}_2,{ℛ}_3)$	242.51	491.02 $\Delta >>10$	498.52 $\Delta >>10$	229.44	464.89 $\Delta >>10$	472.07 $\Delta >>10$
$({ℛ}_1,{ℛ}_3,{ℛ}_4)$	170.58	347.16 $w=0.07$ $\Delta =5.0$	354.66 $w=0.15$ $\Delta =3.0$	144.58	295.15 $w=0.64$	302.34 $w=0.81$
${\left({ℛ}_i\right)}_{i\in ⟦1,4⟧}$	167.05	342.10 $w=0.90$	351.68 $w=0.68$	144.24	296.48 $w=0.33$ $\Delta =1.3$	306.06 $w=0.12$ $\Delta =3.7$

Table 2. Table 2: Parameter sets used to generate the in silico datasets. We use two distinct parameter sets for each submodel, shown in the two rows associated to each submodel.

	${\alpha }_1$	$\beta$	${\alpha }_2$	$\gamma$	$\mu$
$({ℛ}_1,{ℛ}_3)$	1	0	0.7	0	5
	1	0	0.007	0	5
$({ℛ}_1,{ℛ}_4)$	1	0	0	0.7	5
	1	0	0	0.007	5
$({ℛ}_1,{ℛ}_2,{ℛ}_3)$	1	0.01	0.07	0	5
	1	100	0.07	0	5
$({ℛ}_1,{ℛ}_3,{ℛ}_4)$	1	0	0.007	0.7	5
	1	0	0.007	0.07	5
$({ℛ}_1,{ℛ}_2,{ℛ}_4)$	1	0.01	0	0.07	5
	1	100	0	0.07	5

Table 3. Table 3: Size-step used for each parameter in the PLE estimate, in log-scale, within each submodel and each datasets.

Model	Parameter	Experimental Wild-Type/Mutant datasets	in silico Dataset $1$	in silico Dataset 2
$({ℛ}_1,{ℛ}_4)$	$\mu$	0.015	0.005	0.01
	$\gamma$	0.12	0.01	0.06
$({ℛ}_1,{ℛ}_3)$	$\mu$	0.015	0.005	0.005
	${\alpha }_2$	0.04	0.01	0.01
$({ℛ}_1,{ℛ}_2,{ℛ}_4)$	$\mu$	0.015	0.01	0.015
	$\beta$	0.12	0.07	0.12
	$\gamma$	0.12	0.07	0.12
$({ℛ}_1,{ℛ}_2,{ℛ}_3)$	$\mu$	0.015	0.015	0.015
	$\beta$	0.12	0.12	0.12
	${\alpha }_2$	0.12	0.02	0.02
$({ℛ}_1,{ℛ}_3,{ℛ}_4)$	$\mu$	0.01	0.01	0.01
	${\alpha }_2$	0.08	0.01	0.01
	$\gamma$	0.08	0.01	0.01
$({ℛ}_1,{ℛ}_2,{ℛ}_3,{ℛ}_4)$	$\mu$	0.015	0.015	0.015
	$\beta$	0.12	0.12	0.12
	${\alpha }_2$	0.12	0.12	0.12
	$\gamma$	0.12	0.12	0.12

Table 4. Table 4: Wild-Type MLE parameter sets. MLE estimates and confidence intervals for each submodel using the likelihood given by Eqs.( 15 )-( 18 ), and (last row) for the initial condition parameter using likelihood given by Eq. ( 19 ).

Model	$\beta$	${\alpha }_2$	$\gamma$	$\mu$
$({ℛ}_1,{ℛ}_4)$	/	/	${10}^{-6}$ $\in (0;0.12]$	$7.49$ $\in [7.05;7.83]$
$({ℛ}_1,{ℛ}_3)$	/	$1.18$ $\in [0.67;1.57]$	/	$7.22$ $\in [6.81;7.83]$
$({ℛ}_1,{ℛ}_2,{ℛ}_4)$	${10}^6\in ℝ$	/	${10}^{4.35}$ $\in (0;{10}^{5.03}]$	$7.45$ $\in [7.05;7.83]$
$({ℛ}_1,{ℛ}_2,{ℛ}_3)$	${10}^6$ $\in [1.52;+\infty )$	${10}^{5.75}$ $\in [2.00;{10}^{5.88}]$	/	$7.07$ $\in [5.15;6.35]$
$({ℛ}_1,{ℛ}_3,{ℛ}_4)$	/	$0.27$ $\in [0.022;0.52]$	${10}^{-}6$ $\in (0;0.068]$	$7.20$ $\in [6.69;7.69]$
$({ℛ}_1,{ℛ}_2,{ℛ}_3,{ℛ}_4)$	${10}^6$ $\in [4.64;+\infty )$	${10}^{4.78}$ $\in [0.87;{10}^{5.27}]$	${10}^{-6}$ $\in (0;{10}^{4.67}]$	$7.06$ $\in [6.58;7.56]$
Primordial follicle dataset	/	/	/	$6.22$ $\in [5.54;6.67]$

Table 5. Table 5: Mutant parameter sets. MLE estimates and confidence intervals for each submodel using the likelihood given by Eqs.( 15 )-( 18 ), and (last row) for the initial condition parameter using likelihood given by Eq. ( 19 ).

Model	$\beta$	${\alpha }_2$	$\gamma$	$\mu$
$({ℛ}_1,{ℛ}_4)$	/	/	0.14 $\in (0;0.28]$	6.40 $\in [5.93;6.81]$
$({ℛ}_1,{ℛ}_3)$	/	1.63 $\in [1.26;2.20]$	/	5.91 $\in [5.34;6.35]$
$({ℛ}_1,{ℛ}_2,{ℛ}_4)$	${10}^6$ $\in ℝ$	/	${10}^{5.11}$ $\in [0.12;{10}^{5.39}]$	6.26 $\in [5.72;6.81]$
$({ℛ}_1,{ℛ}_2,{ℛ}_3)$	${10}^6$ $\in ℝ$	${10}^6$ $\in [1.52;+\infty )$	/	5.57 $\in [5.15;6.35]$
$({ℛ}_1,{ℛ}_3,{ℛ}_4)$	/	$0.52$ $\in [0.21;0.91]$	${10}^{-6}$ $\in (0,0.98]$	5.94 $\in [5.43;6.54]$
$({ℛ}_1,{ℛ}_2,{ℛ}_3,{ℛ}_4)$	$2.81$ $\in {ℝ}_{+}$	$1.16$ $\in [0.28;{10}^{5.51}]$	${10}^{-6}$ $\in (0;{10}^{4.9}]$	$5.83$ $\in [5.15;6.35]$
Primordial follicle dataset	/	/	/	$6.77$ $\in [5.75;7.60]$

Equations

$$\begin{array}[]{lcl}&\text{ Cell events }&\text{ Rate }\\ \mathcal{R}_{1}:&(F,C)\rightarrow(F-1,C+1),&\,\,\,\alpha_{1}F,\\ \mathcal{R}_{2}:&(F,C)\rightarrow(F-1,C+1),&\,\,\,\beta\frac{FC}{F+C},\\ \mathcal{R}_{3}:&(F,C)\rightarrow(F,C+1),&\,\,\,\alpha_{2}F,\\ \mathcal{R}_{4}:&(F,C)\rightarrow(F,C+1),&\,\,\,\gamma C\,.\end{array}$$

$$\begin{gathered} F_t=F_0-{\mathcal{Y}}_1\left({\int }_0^t{\alpha }_1{F}_{s}ds\right)-{\mathcal{Y}}_2\left({\int }_0^t\beta \frac{F_s{C}_s}{F_s+C_s}ds\right), \\ {C}_t={\mathcal{Y}}_1\left({\int }_0^t{\alpha }_1{F}_{s}ds\right)+{\mathcal{Y}}_2\left({\int }_0^t\beta \frac{F_s{C}_s}{F_s+C_s}ds\right) \\ +{\mathcal{Y}}_3\left({\int }_0^t{\alpha }_2{F}_{s}ds\right)+{\mathcal{Y}}_4\left({\int }_0^t\gamma C_{s}ds\right). \end{gathered}$$

$$\begin{gathered} \mathcal{L}g(f,c)=({\alpha }_{1}f+\beta \frac{fc}{f+c})\left[g(f-1,c+1)-g(f,c)\right] \\ +({\alpha }_{2}f+\gamma c)\left[g(f,c+1)-g(f,c)\right], \end{gathered}$$

$$M_t^g:=g(X_t)-g(X_0)-{\int }_0^t\mathcal{L}g(X_s)ds$$

$$\tau :=\inf \{t\ge 0;F_t=0|F_0\}.$$

$${\tau }_L:=\inf \{t;F_t^L=0|f_0\}.$$

$$T_{k+1}:=T_k+\mathcal{E}\left({\alpha }_1(f_0-k)\right),$$

$$f_{{\tau }_L}(t)={\alpha }_1{f}_0{e}^{-{\alpha }_{1}t}(1-e^{-{\alpha }_{1}t}{)}^{f_0-1}\text{mathds}{1}_{[0,+\infty [}(t),$$

$$C_t^{0,0}=1+\mathcal{Y}\left(\gamma {\int }_0^t{C}_s^{0,0}ds\right),$$

$$C_t^L=\underset{\text{ cell lineages generated by cell event }{\mathcal{R}}_1}{\underbracket{\sum _{k=1}^{F_0}{C}_{t-T_k^0}^{k,0}\text{mathds}{1}_{t\ge T_k^0}}}+\underset{\text{cell lineages generated by cell event }{\mathcal{R}}_3}{\underbracket{\sum _{k=0}^{F_0-1}\sum _{j=1}^{N_k(t)}{C}_{t-T_k^j}^{k,j}\text{mathds}{1}_{t\ge T_k^j}}},$$

$$N_k(t)={\mathcal{Y}}_3\left({\alpha }_2{\int }_0^{t\wedge T_{k+1}}{F}_s^{L}ds\right)-{\mathcal{Y}}_3\left({\alpha }_2{\int }_0^{T_k}{F}_s^{L}ds\right).$$

$$T_k^j:=T_k^{j-1}+\mathcal{E}\left({\alpha }_2(f_0-k)\right),$$

$$N_k(t)=\sum _j{1}_{t<\min (T_k^j,T_{k+1})}.$$

$$\mathbb{E}\left[C_{{\tau }_L}^L\right]<\infty$$

$$\mathbb{E}\left[C_{{\tau }_L}^L\right]=f_0(1+\frac{{\alpha }_2}{{\alpha }_1}).$$

$$\mathbb{E}\left[C_{{\tau }_L}^L\right]=1+\sum _{k=1}^{f_0-1}\frac{(f_0-k)!}{\left((f_0-k)-\frac{\gamma }{{\alpha }_1}\right)!}$$

$$\mathbb{E}\left[{\left(C_{\tau }\right)}^p\right]<\infty ,$$

$$p\gamma <{\alpha }_1+\beta .$$

$${\tau }_n=\min \left(\inf \{t\ge 0,C_t\ge n\},\tau \right),$$

$${\mathcal{L}}_F^n\varphi (f)=\left({\alpha }_1+\beta \frac{n}{1+n}\right)f\left[\varphi (f-1)-\varphi (f)\right],$$

$${\mathcal{L}}_C\varphi (c)=\left[({\alpha }_1+\beta +{\alpha }_2)f_0+\gamma c\right]\left[\varphi (c+1)-\varphi (c)\right],$$

$$F_t^n=f_0,C_t^n=n$$

$$C_{\tau }\le C_{{\tau }_n+{\overline{\tau }}_n}^n.$$

$${\overline{\tau }}_n:=\inf \{t>0;F_{t+{\tau }_n}^n=0|f_0\}.$$

$$F_t\le F_t^n,C_t\le C_t^n.$$

$$F_{t+{\tau }_n}\le F_{t+{\tau }_n}^n,C_{t+{\tau }_n}\le C_{t+{\tau }_n}^n,\text{a.s.}$$

$$F_{{\tau }_n}\le F_{{\tau }_n}^n,C_{{\tau }_n}\le C_{{\tau }_n}^n.$$

$$\left({\alpha }_1+\beta \frac{C_{t+{\tau }_n}}{F_{t+{\tau }_n}+C_{t+{\tau }_n}}\right)F_{t+{\tau }_n}\ge \left({\alpha }_1+\beta \frac{n}{1+n}\right)F_{t+{\tau }_n},$$

$$({\alpha }_1+{\alpha }_2)F_{t+{\tau }_n}+\beta \frac{C_{t+{\tau }_n}}{F_{t+{\tau }_n}+C_{t+{\tau }_n}}{F}_{t+{\tau }_n}\le ({\alpha }_1+{\alpha }_2+\beta )f_0,$$

$$C_{\tau }=C_{\tau -{\tau }_n+{\tau }_n}\le C_{\tau -{\tau }_n+{\tau }_n}^n\le C_{{\overline{\tau }}_n+{\tau }_n}^n,$$

$$p\gamma <{\alpha }_1+\beta \frac{n}{1+n}$$

Full text

1 Introduction

In mammals, the number of oocytes (egg cells) available for a female throughout her reproductive life is fixed once for all, during the fetal or perinatal period monniaux_18 . Dormant oocytes are sheltered within somatic structures called ovarian follicles, which remain in a quiescent state until they get activated and undergo a longstanding process of growth and maturation ending by ovulation (release of a fertilizable oocyte). Growth initiation is asynchronous among follicles, so that all developmental stages can be observed in the ovaries at a given time, and follicles can remain quiescent for as long as tens of years reddy_10 .

In the earliest stages of development, ovarian follicles are made up of the oocyte and a single layer of surrounding somatic cells. The initial cell number is on the order of ten or several of tens according to the species and is quite variable between follicles. Such a variability is inherited from the mechanism underlying the formation of primordial follicles monniaux_18b ; sawyer_02 , which assemble from the fragmentation of syncytium structures (the germ cell cysts) and retrieve more or less somatic cells.

The activation of primordial (quiescent) follicles is characterized by three main processes picton_01 : (i) an irreversible transition of the somatic cell phenotype, characterized by a change in their shape, from flattened (precursor cells) to cuboidal (proliferative cells); (ii) an increase in the number of somatic cells by cell division and (iii) the awakening and associated enlargement of the oocyte. The activation phase is ended when all somatic cells have transitioned, at which time the mono-layer developmental stage is completed, and somatic cells will go on proliferating and build up several concentric layers fortune_03 ; clement_coupled_2013 ; CRY2019 .

In this work, we focus on the sequence of events occurring just after the initiation of follicle growth. A key issue is to determine whether cell proliferation is concomitant or posterior to cell shape change, and to assess both the time needed for all precursor cells to complete transition and the corresponding increase in the cell number with respect to the initial cell number.

We introduce a continuous-time Markov chain model for cell population dynamics accounting for both cell transition and division. Within such a formalism, linear models have been built up on the branching property, disregarding cellular interactions kimmel_theory_1963 ; harris_theory_1963 , while nonlinear models have accounted for interactions among different cell populations (e.g., typically, a feedback from differentiated cells onto precursor cells) either to ensure homeostasis, as in dynamical models for blood cells getto_mathematical_2015 ; stiehl_stem_2017 ; pujo_blood_2016 , or to achieve a proper developmental sequence, as in dynamical models for neural cells freret-hodara_16 . On our side, we are interested in assessing the duration of the activation process, i.e. the extinction time of the population of precursor cells, and in ordering the events taking place during activation. A natural concept in probability theory to investigate these issues is the first passage time theory darling_first_1953 ; van_kampen_stochastic_1992 , which aims to characterize the statistics of random events related to some particular outcomes. The analysis of first passage times are becoming more and more popular in mathematical biology chou_first_2014 ; castro_mathematical_2015 , to quantify random times needed to reach a given final state, such as population extinction for instance.

Typically, the parameters of cell dynamics models are calibrated using time series of cell counts sorted into different cell types marr_multiscale_2012 ; glauche_lineage_2007 . In contrast, in the case of early folliculogenesis, precursor and proliferative cell numbers are not available directly as a function of time, but only in relation with other morphological variables such as the oocyte and follicle diameters braw-tal_studies_1997 ; gougeon_morphometric_1987 ; lundy_populations_1999 ; meredith_classification_2000 , so that we lack kinetic information. To overcome this difficulty, we use the embedded discrete-time Markov chain to apply classical statistical tools like the maximum likelihood wilkinson_markov_2013 , and parameter identifiability concepts raue_structural_2009 .

The manuscript is organized as follows. In Section 2, we introduce a stochastic model of cell population dynamics, with two state variables and four cell events (reactions). In section 3 we analyze both the linear and nonlinear versions of the model in the Markov chain framework. In the linear case, we obtain analytical formulas for the mean extinction time. In the nonlinear case, we design a numerical scheme based on a rigorous Finite State Projection (see munsky_finite_2006 ; kuntz_deterministic_2017 ) and coupling techniques to assess the mean extinction time. In both cases, we study the sensitivity of the extinction time, as well as of the proliferative cell number at extinction time, with respect to the parameter values. In section 4, using the embedded Markov chain, we calibrate the parameters of the model from experimental, time-free datasets, and analyze the practical identifiability. Using model selection criteria, we identify the key parameters that shed light on the most likely events that occur during the activation process. From this data-fitting approach, we manage to retrieve hidden kinetic information and provide some biological interpretations of our results. We conclude in section 5.

2 Model design and formulation

Our model allows us to study the joint dynamics of the precursor cells $F$ and proliferative cells $C$ within a single follicle, whose populations are ruled by four types of possible cell events. In the absence of specific information, we used the simplest formulation as possible for all event rates, according to Occam’s razor principle. Two cell events occur at the expense of the precursor cells, which are consumed during their transition : (i) $\mathcal{R}_{1}$ is the spontaneous transition of precursor cells into proliferative cells, whose rate $\alpha_{1}F$ is linearly proportional to the number of precursor cells; (ii) $\mathcal{R}_{2}$ is the auto-amplified transition of precursor cells into proliferative cells, which occurs at rate $\beta\frac{FC}{F+C}$. This event represents the feedback of proliferative cells onto the transition of the precursor cells. Two other cell events increase the proliferative cell population without affecting the precursor cell population: (i) $\mathcal{R}_{3}$ is an asymmetric division of precursor cells $F$ (giving rise to one precursor cell and one proliferative cell), which occurs at rate $\alpha_{2}F$; (ii) $\mathcal{R}_{4}$ is a symmetric division of the proliferative cells $C$ (giving rise to two proliferative cells), which occurs at rate $\gamma C$.

These four cell events are the building blocks of the main model $\mathcal{M}_{FC}$, which is summarized below :

[TABLE]

Cell events $\mathcal{R}_{1}$ and $\mathcal{R}_{4}$ constitute the fundamental ingredients involved in the activation process. We also consider two additional cell events, $\mathcal{R}_{3}$ and $\mathcal{R}_{{2}}$, which are not only intended to enrich the model behavior, but are also substantiated by biological observations.

Cell event $\mathcal{R}_{1}$ corresponds to the spontaneous transition undergone by a precursor cell, including the very first event. Firing only $\mathcal{R}_{1}$ events is sufficient to complete activation, yet in this case the final cell number is unchanged with respect to the initial number, which is not what is systematically observed in the experimental data lundy_populations_1999 ; gougeon_morphometric_1987 ; lintern_79 . On the scale of a whole follicle, the awakening of precursor cells triggers the exit from the primordial follicle pool and initiate the process of follicle growth and development zhang_somatic_2014 . Awakening is induced by activation of the protein complex mTORC1 in somatic cells (and not oocyte), by oxygen and stress-or energy-induced metabolites reaching the follicle environment in the ovarian cortex zhang_somatic_2014 . A natural choice for this spontaneous reaction is to consider that it occurs independently in each precursor cell, so that the transition rate of precursor cells is proportional (with coefficient $\alpha_{1}$) to their number $F$.

Cell event $\mathcal{R}_{2}$ corresponds to auto-amplified precursor cell transitions. It has the same cell output as cell event $\mathcal{R}_{1}$ (loss of one precursor cell), yet it can speed up the transition rate after the first triggering event. The amplification is mediated by a positive feedback exerted by already transitioned, proliferative cells; the transition rate is [math] when $C=0$, and increases with $C$. Such an auto-amplification is expected to result from the molecular mechanisms underlying follicle activation and establishing a dialog between the oocyte and somatic cells monniaux_16 . The activated somatic cell(s) start stimulating the oocyte through specific signaling pathways (KIT-Ligand cytokine). In turn, once activated, the oocyte signals to the somatic cells through factors of the TGF$\beta$ family knight_06 (mainly GDF9 and BMP15). This molecular dialog settles a positive feedback loop, which can be represented by an auto-amplified transition rate. In sheep, there exist natural mutations affecting this molecular dialog (disruption of either the GDF9 or BMP15 ligand, or the receptor to BMP15). Introducing cell event $\mathcal{R}_{2}$ can help to investigate possible differences in the activation process in wild-type compared to mutant strains. More specifically, we have access to experimental cell numbers (courtesy of Ken McNatty) obtained either from a wild-type strain (Ile-de-France) or a mutant strain for BMP15R (Booroola), whose follicle development is known to be clearly different in the multi-layer stages lundy_populations_1999 , especially as far as cell dynamics. Whether cell dynamics is also affected during the mono-layer stage remains unclear reader_12 , which is an additional motivation for this work. In the analysis performed in Section 3, the specific formulation of the reaction rate of event $\mathcal{R}_{2}$ does not matter much. We just need to assume that it is linearly bounded by $F$, which is sensible with respect to cell cycle constraints. To fit the model to available data in Section 4, we needed to specify further the shape of the nonlinearity, and chose a parameterization including as few parameters as possible. We refer to Appendix 6.1 for a basic justification of this choice.

Cell event $\mathcal{R}_{3}$ corresponds to a self-renewal transition event that does not consume a precursor cell and may be fired as a first triggering event. It is directly inspired from asymmetric division events commonly observed in developmental cell lineages. The speculation that precursor (flattened) cells might divide while transitioning is compatible with experimental studies where KI67 staining (a marker of cell cycle progression) was detected in some flattened cells dasilva_08 . Since the number of flattened cells is non increasing, one can envisage the existence of self-renewing asymmetric divisions in flattened cells, giving birth to one proliferative cell (and keeping the precursor cell number unchanged). The rate of event $\mathcal{R}_{3}$ is chosen in a similar way as that of event $\mathcal{R}_{1}$.

Cell event $\mathcal{R}_{4}$ corresponds to the symmetric division of proliferative cells, giving birth to two identical daughter cells. It has the same cell output as cell event $\mathcal{R}_{4}$ (gain of one proliferative cell). All $C$ cells are supposed to progress through the cell cycle (growth fraction of one), and there is no cell loss at mitosis, so that the proliferative population growths exponentially (with rate $\gamma$).

All the reactions rates ($\alpha_{1}$, $\beta$, $\alpha_{2}$ and $\gamma$) are non-negative. At initial time, there are only precursor cells, and the initial condition is chosen as a random positive integer variable, in consistency with the observed biological variability.

In the following, we will use different submodels derived from the full model $\mathcal{M}_{FC}$, by removing either one or several cell events (hence setting to zero the corresponding parameter values $\beta$, $\alpha_{2}$ and/or $\gamma$). We will name these submodels by explicitly mentioning the remaining events. For instance, model ($\mathcal{R}_{1},\mathcal{R}_{3}$) consists only of the spontaneous cell transition event and asymmetric cell division ($\beta=\gamma=0$), while model ($\mathcal{R}_{1},\mathcal{R}_{4}$) is composed of the spontaneous cell transition event and asymmetric cell division ($\beta=\alpha_{2}=0$).

Model $\mathcal{M}_{FC}$ is mathematically formulated as a Continuous time Markov chain (CTMC). The stochastic description is especially appropriate when dealing with a small number of cells. Before introducing a precise mathematical formulation, we can illustrate the dynamics of both the precursor and proliferative cells (Figure 1). Initially, the whole population is made of precursor cells ($C=0$), so that the first event has to be a triggering event generating the first proliferative cell ($\mathcal{R}_{1}$ or $\mathcal{R}_{3}$). The $C$ population grows as the $F$ population decreases until extinction (top-left panel). The last event before extinction of precursor cells has to be a consuming event ($\mathcal{R}_{1}$ or $\mathcal{R}_{2}$). The proportion of proliferative cells $p_{C}:=\frac{C}{F+C}$ increases monotonously from [math] to $1$ (bottom-left panel). In the $(C,F)$ phase plane (top-right panel), we can observe that the number of precursor cells remains constant (aligned red or black points on the horizontal line $(k,F)$, $k\in\mathbb{N}$) whenever there is a division event ($\mathcal{R}_{3}$ or $\mathcal{R}_{4}$). In contrast, whenever there is a transition event ($\mathcal{R}_{1}$ or $\mathcal{R}_{2}$), the number of precursor cells decreases by one, as illustrated by the jump from the current line ($(k,F),k\in\mathbb{N}$) to the lower one ($(k,F-1)$, $k\in\mathbb{N}$). Hence, in this simulation, we observe a sequence of transition and division events (which appear to be here mainly spontaneous transitions $\mathcal{R}_{1}$ and asymmetric divisions $\mathcal{R}_{3}$ due to the specific parameter choice). If we are only given the sequence of events in this plane, we cannot discriminate $\mathcal{R}_{1}$ from $\mathcal{R}_{2}$, neither $\mathcal{R}_{3}$ from $\mathcal{R}_{4}$. Note that, depending on the initial condition, some parts of the phase plane cannot be reached. The trajectories can also be observed in the $(C,p_{C})$ phase plane (bottom-right panel). In this case, the trajectories remain on the curves parameterized by ($(k,\frac{k}{F+k})$, $k\in\mathbb{N}$) if a division event ($\mathcal{R}_{3}$ or $\mathcal{R}_{4}$) occurs, whereas they move to the upper curves parameterized by ($(k,\frac{k}{F-1+k})$, $k\in\mathbb{N}$) whenever a transition event ($\mathcal{R}_{1}$ or $\mathcal{R}_{2}$) occurs.

Model formulation and hypotheses

On a probability space $(\Omega,\mathcal{F},\mathbb{P})$, let the initial number of flattened cells $F_{0}$ be a positive integer random variable. The population of precursor cells $F$ and proliferative cells $C$ follows the Stochastic Differential Equation (SDE) below:

[TABLE]

where $\mathcal{Y}_{i}$, for all $i=1,2,3,4$, are mutually independent standard Poisson processes. $X=(X_{t})_{t\geq 0}$, with $X_{t}:=(F_{t},C_{t})$ for all $t\geq 0$, denotes the solution of (1). $(\mathcal{F}_{t})_{t\geq 0}$ denotes the canonical filtration generated by the process X.

Classically, X can also be seen as a continuous-time Markov chain with countable state space $\mathcal{S}:=\mathbb{N}^{2}\backslash\{(0,0)\}$ whose infinitesimal generator $\mathcal{L}$ is given by

[TABLE]

for all $g$ bounded functions and for all $(f,c)\in\mathcal{S}$.

In the whole study, we will need the following hypotheses:

Hypothesis 1

The spontaneous activation rate $\alpha_{1}$ is positive.

Hypothesis 2

The initial condition $F_{0}$ is $L_{2}$-integrable.

With Hypothesis 2, we apply Theorem 1.22 of anderson_stochastic_2015 (p.12-13) and deduce that the process $M^{g}_{t}$ defined as

[TABLE]

is a $\mathcal{F}_{t}$-martingale, for all $t\geq 0$ and any bounded function $g$.

Note that process $F$ is a non-negative decreasing process. To study the hitting time of the state $F=0$, we introduce the following definition

Definition 1

Let $\tau$ be the extinction time of the precursor cell population $F$

[TABLE]

The number of proliferative cells $C$ at $t=\tau$ is $C_{\tau}$.

To control the first moment of $C_{\tau}$, the number of proliferative cells at the extinction time, we will also need an additional hypothesis:

Hypothesis 3

The maximal activation rate $\alpha_{1}{+\beta}$ is strictly greater than the proliferation rate $\gamma$: $\alpha_{1}{+\beta}>\gamma$.

3 Model analysis

In this section we analyze the mean extinction time of the precursor cell population and the number of proliferative cells at extinction. We start in subsection 3.1 by recalling some analytical formulas for model ($\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4}$) (when $\beta=0$, linear rate functions). Then, in subsection 3.2, we deduce a necessary and sufficient condition to ensure that the mean of $C_{\tau}$ is finite for the complete model ($\mathcal{M}_{FC}$), by finding a lower and an upper-bound thanks to a coupling argument. In subsection 3.3, we finally use the upper-bound in a finite-state projection algorithm to obtain an efficient way to simulate the means of $\tau$ and $C_{\tau}$, and numerically investigate the role of the feedback rate $\beta$ on their values.

To simplify the proofs, we will consider in the following that the initial condition is a deterministic value $f_{0}\in\mathbb{N}^{*}$. All the proofs can be generalized to the random $F_{0}$ case by conditioning by the law of $F_{0}$.

3.1 Analytical expressions in the linear case, model ($\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4}$).

When $\beta=0$, process $X$ is linear, and we can compute the law of the extinction time. In the case when, in addition, $\alpha_{2}=0$ and/or $\gamma=0$, the mean of $C_{\tau}$ can also be computed.

In this subsection we will write $X^{L}_{t}=(F^{L}_{t},C^{L}_{t})$ the solution of SDE (1) when $\beta=0$ and $\tau_{L}$ the associated extinction time of the population $F^{L}_{t}$:

[TABLE]

Note that process $F^{L}$ is independent of process $C^{L}$. The jumping times $T_{k}$ of $F^{L}$, for all $k\in\llbracket 0,f_{0}-1\rrbracket$, are given by

[TABLE]

with $T_{0}=0$ by convention, and $\mathcal{E}(\lambda)$ denotes an exponential random variable of mean $1/\lambda$. Note that $T_{f_{0}}=\tau_{L}$.

Proposition 1 ($F^{L}_{t}$ and $\tau_{L}$ laws)

Under Hypothesis 1 and for all $t\geq 0$, $F^{L}_{t}$ follows a binomial law with parameters $(f_{0},e^{-\alpha_{1}t})$, and the extinction time $\tau_{L}$, defined by formula (3), follows a generalized Erlang law (or hypo-exponential law) of density:

[TABLE]

such that $\mathbb{E}\left[\tau_{L}\right]=\displaystyle\frac{1}{\alpha_{1}}\sum_{k=1}\frac{1}{k}$.

The proof of Proposition 1 is standard and given in Appendix 6.3 for the reader convenience.

We now study the mean number of proliferative cells at the extinction time. We first decompose process $C_{t}^{L}$ as a sum of elementary processes. We introduce the following binary branching process

[TABLE]

where $\mathcal{Y}$ is a Poisson process. This process is often referred as the Yule process (Bailey, , Chap. 8). We then define the stochastic processes $C^{k,j}$, for $(k,j)\in\mathbb{N}\times\mathbb{N}$, as independent and identically distributed Yule processes. Each process $C^{k,j}$ represents the cell population that arises from the successive symmetric divisions starting from a single newly transitioned proliferative cell. We thus refer the $C^{k,j}$ processes as ”cell lineages”. Process $C^{L}_{t}$ is a branching process with immigration driven by cell events $\mathcal{R}_{1}$ and $\mathcal{R}_{3}$, it can indeed be written as the sum of the cell lineages $C^{k,j}$ (illustrated in Figure 2, such lineage decomposition goes back at least to clonal cell population size studies like Luria_43 ): for all $t\geq 0$,

[TABLE]

where we define, for all $k\in\llbracket 1,f_{0}\rrbracket$,

• •

$T^{0}_{k}:=T_{k}$ (with $T_{k}$ given by equation (4)), the $k$-th jumping time of cell event $\mathcal{R}_{1}$.

• •

$N_{k}(t)$, the number of occurrences of cell event $\mathcal{R}_{3}$ between $T_{k}$ and $T_{k+1}$, for $t\geq T_{k}$. Note that

[TABLE]

• •

for all $j\in\llbracket 1,N_{k}(t)\rrbracket$,

[TABLE]

the $j$-th jumping time of the cell event $\mathcal{R}_{3}$ occurring between random times $T_{k}$ and $T_{k+1}$. Hence

[TABLE]

We use the sum defined in Eq. (6) to obtain a necessary and sufficient condition to ensure that the mean of $C^{L}_{\tau_{L}}$ is finite, and to obtain analytical formulas for the submodels $(\mathcal{R}_{1},\mathcal{R}_{3})$ ($\beta=\gamma=0$) and $(\mathcal{R}_{1},\mathcal{R}_{4})$ ($\beta=\alpha_{2}=0$).

Proposition 2 (First moment of $C^{L}_{\tau_{L}}$)

Under Hypothesis 1 we have that

[TABLE]

if, and only if, Hypothesis 3 holds. Moreover,

For the submodel $(\mathcal{R}_{1},\mathcal{R}_{3})$ ($\gamma=0$), we have

[TABLE] 2. 2.

For the submodel $(\mathcal{R}_{1},\mathcal{R}_{4})$ ($\alpha_{2}=0$), we have

[TABLE]

where we used notation $(m-x)!=\prod_{i=1}^{m}(i-x)$ for $m\in\mathbb{N}^{*}$ and $x\in[0,1)$.

The proof of Proposition 2 is classical and is given in Appendix 6.3 for the reader convenience.

Remark 1

A simple analytical formula cannot be obtained for the first moment of $C^{L}_{\tau_{L}}$ for submodel ($\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4}$) since it is tricky to deal with expectation in the second term of relation (6).

3.2 Lower and upper bounds of the nonlinear model ($\mathcal{M}_{FC}$)

In the general case, we cannot obtain analytical expressions for $\tau$, and we will rather use numerical simulations. To control the numerical error, we need tractable bounds of the stochastic model introduced in Eq. (1), which are obtained in this subsection. We first note that all moments of $\tau$ are unconditionally finite, as $F_{t}$ decreases by one at rate at least $\alpha_{1}F_{t}$, so that $\tau$ is stochastically dominated by $\tau_{L}$ given in Proposition 1. Our main result is a necessary and sufficient condition to obtain finite moments for $C_{\tau}$.

Theorem 3.1

For any $p\geq 1$, we have, under Hypothesis 1, that

[TABLE]

if, and only if

[TABLE]

In particular, for $p=1$, $\mathbb{E}\left[\left(C_{\tau}\right)\right]<\infty$ if and only if Hypothesis 3 holds.

The main step in the proof of Theorem 3.1 is Proposition 3, which provides us with a concrete upper bound usable to control the numerical error of our algorithm. We first need the following definitions to set up the upper-bound process:

Definition 2 (truncated extinction time)

Let $X_{t}:=(F_{t},C_{t})$, for all $t\geq 0$, the solution of (1). For any $n\geq 1$, we define

[TABLE]

Definition 3 (upper-bound process)

For any $n\geq 1$, we define the following infinitesimal generators

[TABLE]

for any $\phi$ bounded on $\mathbb{N}$, and any $f,c\in\mathbb{N}$.

Proposition 3 (Upper-bound)

Let $X_{t}:=(F_{t},C_{t})$, for all $t\geq 0$, the solution of (1). For any $n\geq 1$, there exists a couple $(F^{n},C^{n})$ such that, for any $t<\tau_{n}$

[TABLE]

*and $(F_{t+\tau_{n}}^{n},C_{t+\tau_{n}}^{n})_{t\geq 0}$ is a continuous time Markov chain of generators

$(\mathcal{L}^{n}_{F},\mathcal{L}_{C})$ satisfying,*

[TABLE]

where

[TABLE]

The proof of Proposition 3 proceeds by a coupling argument between $(F,C)$ and $(F^{n},C^{n})$. The random variable $C_{\tau_{n}+\overline{\tau}_{n}}^{n}$ has a finite $p$-moment under (9), and this moment is analytically tractable, thanks to Proposition 6 in Appendix 6.3.

Proof (of Proposition 3)

Let $F_{t},C_{t}$ given by equation (1). Let $n\geq 1$ and $(F^{n},C^{n})$ defined on $(0,\tau_{n})$ by (10). Clearly, for $t\in(0,\tau_{n})$,

[TABLE]

Then, we may choose $(F_{t+\tau_{n}}^{n},C_{t+\tau_{n}}^{n})_{t\in\mathbb{R}^{+}}$ such that its infinitesimal generator is given by $(\mathcal{L}^{n}_{F},\mathcal{L}_{C})$ in Definition 3, and such that its trajectories satisfy, for any $t\geq 0$,

[TABLE]

This is clearly possible as

For $t=0$, the initial condition satisfies

[TABLE]

For any $t>0$, we can that ensure $F_{t+\tau_{n}}$ stays below $F_{t+\tau_{n}}^{n}$ because

[TABLE]

as $C_{t}$ is non-decreasing.

For any $t>0$, we can ensure that $C_{t+\tau_{n}}$ stays below $C_{t+\tau_{n}}^{n}$ because

[TABLE]

as $F_{t}$ is non-increasing.

On the event $(\tau_{n}<\tau)$, we have $C_{t+\tau_{n}}\leq C_{t+\tau_{n}}^{n}$ for all times $t\geq 0$, and $0<\tau-\tau_{n}\leq\overline{\tau}_{n}$. We note that $\overline{\tau}_{n}$ follows a generalized Erlang law of parameter $\alpha_{1}+\beta\frac{n}{1+n}$, by straightforward adaptation of Proposition 1. Hence

[TABLE]

as $C^{n}(t)$ is non-decreasing.

On the event $(\tau_{n}=\tau)$, we clearly have $C_{\tau}\leq n\leq C_{\overline{\tau}_{n}+\tau_{n}}^{n}$.

Combining both cases proves equation (11).

$\square$

We now proceed to the proof of Theorem 3.1.

Proof (of Theorem 3.1)

Let $p\geq 1$. Condition (9) implies that there exists (a sufficiently large) $n\geq 1$ such that

[TABLE]

holds true. In particular, for such a $n$, we have

[TABLE]

by consequence of Proposition 6 in Appendix 6.3 (which also provides explicit bounds for $p=1,2$). We now prove that condition (9) is necessary, excluding the trivial situation where $f_{0}=1$ and $\alpha_{2}=0$ (in which case $C_{\tau}=1$). First, note that $\tau$ is stochastically lower-bounded by an exponential random variable of rate $(\alpha_{1}+\beta)$, because the maximal activation rate of $F$ is $\alpha_{1}+\beta$. The Yule process in Eq. (5) provides a lower-bound for $C_{t}$ for times $t$ greater than the first event (given by an exponential random variable of rate $\alpha_{1}+\alpha_{2}$). Thus, the Yule process stopped at an exponential time of parameter $\alpha_{1}+\beta$ provides a lower bound for $C_{\tau}$. We conclude again by Proposition 6 in Appendix 6.3 (with $n=1$).

$\square$

3.3 Numerical scheme for the mean extinction time and mean number of proliferative cells at the extinction time

We now have all the ingredients to study numerically the impact of the model parameters on the mean activation duration of an ovarian follicle (mean extinction time of precursor cells) and the mean number of proliferative cells produced during this phase.

From the martingale problem (2), it is a standard result to compute the moment of $\tau$ and $C_{\tau}$. Let the domain $\mathcal{D}$ be defined as

[TABLE]

We look for the value $g(f_{0},0)$ where $g$ is solution of

[TABLE]

where function $g_{0}$ and scalar $\alpha$ are to be chosen according to whether we want to obtain $\mathbb{E}\left[\tau\right]$ or $\mathbb{E}\left[C_{\tau}\right]$.

For $\mathbb{E}\left[\tau\right]$, we take, for all $c\in\mathbb{N}$, $g_{0}(c)=0$ and $\alpha=-1$. 2. 2.

For $\mathbb{E}\left[C_{\tau}\right]$, we take, for all $c\in\mathbb{N}$, $g_{0}(c)=c$ and $\alpha=0$.

We can notice that system (13), which is similar to the Kolmogorov backward equation, is unclosed, and there exists no analytical solution. We now obtain a numerical estimate for the scalar $g(f_{0},0)$ using a domain truncation method, as proposed in munsky_finite_2006 ; kuntz_deterministic_2017 .

Domain truncation method

For $r\in\mathbb{N}^{*}$, let $\mathcal{D}^{r}$ be the following truncated domain

[TABLE]

Note that the truncated extinction time111Although the cut-off $r$ plays a similar role as the index $n$ from section 3.2, we will need two distinct values for the numerical scheme, so that we stick with two different notations, to avoid possible confusion. $\tau_{r}$ defined in Definition 2 is the first exit time from $\mathcal{D}^{r}$,

[TABLE]

As $\mathcal{D}^{r}\subset\mathcal{D}$, we clearly have $\tau_{r}\leq\tau$ and consequently $C_{\tau_{r}}\leq C_{\tau}$. Also, as $\mathcal{D}^{r}$ is a strictly increasing sequence of sets such that $\cup_{r}\mathcal{D}^{r}=\mathcal{D}$, the upper-bound obtained in Proposition 3 will allow us to prove that

[TABLE]

and

[TABLE]

and to control the speed of convergence.

Proposition 4 (Domain truncation relative error)

Let $p\in\mathbb{N}^{*}$, such that $\mathbb{E}[(C_{\tau})^{p}]<\infty$. Then, we have

[TABLE]

and

[TABLE]

where $\overline{\tau}_{r}$ and $C^{r}$ are defined in Proposition 3 (see Eqs. (10)-(12)).

Proof

At $\tau_{r}$, we either have $F_{\tau_{r}}=0$, in which case $\tau_{r}=\tau$, or $F_{\tau_{r}}\in[1,f_{0}]$, in which case $\tau_{r}<\tau$ and $C_{\tau}\geq r$. We use this dichotomy to compute the difference $\mathbb{E}\left[\tau-\tau_{r}\right]$:

[TABLE]

Given that $F_{\tau_{r}}\geq 1$, we have $C_{t}\geq r$ for all $t\in[\tau_{r},\tau]$. Hence, by the same coupling procedure as in the proof of Proposition 3, $\tau-\tau_{r}\leq\overline{\tau}_{r}$ where $\overline{\tau}_{r}$ is independent of $F_{\tau_{r}}$, and follows a generalized Erlang law of parameter $\alpha_{1}+\beta\frac{r}{1+r}$. Moreover, we clearly have

[TABLE]

We then conclude by Chebychev inequality that

[TABLE]

Using the same reasoning we obtain

[TABLE]

and, given that $F_{\tau_{r}}\geq 1$,

[TABLE]

where again $C^{r}$ is defined in Proposition 3 and is independent of $F_{\tau_{r}}$.

$\square$

The sequence of random variables $\overline{\tau}_{r}$ is uniformly bounded in mean,

[TABLE]

Proposition 4 can then be used effectively with any $p\geq 1$ to approximate $\mathbb{E}\left[\tau\right]$ with $\mathbb{E}\left[\tau_{r}\right]$. Using Proposition 3, we deduce that for any $\gamma<\alpha_{1}+\beta$ (e.g. under Hypothesis 3), $\mathbb{E}\left[\tau-\tau_{r}\right]$ is decreasing as $\frac{1}{r}$, with a computable pre-factor given in Proposition 6 in Appendix 6.3.

Similarly, under Hypothesis 3 and from Proposition 6, we obtain the explicit bound for $r$ sufficiently large (such that $\gamma<\alpha_{1}+\beta\frac{r}{1+r}$),

[TABLE]

The latter expression increasing at most linearly in $r$, Proposition 4 can be used with any $p\geq 2$. Under the assumption that $2\gamma<\alpha_{1}+\beta$, we deduce from Proposition 3 that $\mathbb{E}\left[C_{\tau}\right]-\mathbb{E}\left[C_{\tau^{r}}\right]$ is also decreasing as $\frac{1}{r}$, with a computable pre-factor given in Proposition 6 in Appendix 6.3.

To be complete, we detail in Appendix 6.4 the pseudo-code (Algorithm 1) to compute $\mathbb{E}\left[\tau_{r}\right]$ and $\mathbb{E}\left[C_{\tau_{r}}\right]$.

We simulate Algorithm 1 to explore the influence of the amplification rate $\beta$ on both the mean of $\tau$, $\mathbb{E}\left[\tau\right]$, and the mean of $C_{\tau}$, $\mathbb{E}\left[C_{\tau}\right]$ for the nonlinear model.

First, we can prove that both $\mathbb{E}\left[\tau\right]$ and $\mathbb{E}\left[C_{\tau}\right]$ are monotonously decreasing with increasing rate $\beta$, and that the following limits hold, with fixed $(f_{0},\alpha_{1},\alpha_{2},\gamma)$ and under Hypothesis 3:

[TABLE]

The limit $\beta\to 0$ is the consequence of the continuity of the Markov chain with respect to its parameters. For the limit $\beta\to\infty$, note that the very first event to occur is either $\mathcal{R}_{1}$ or $\mathcal{R}_{3}$, and the remaining ones are $\mathcal{R}_{2}$ almost surely in the limit $\beta\to\infty$. These limits are illustrated on Figure 3. Moreover, the different parameter configuration used in Figure 3 leads to the guess222we are not able to prove it, as no analytical formula is available for the full model that both $\mathbb{E}\left[\tau\right]$ (left panel) and $\mathbb{E}\left[C_{\tau}\right]$ (right panel) have a high sensitivity to $\beta$ in the range $\beta\approx\alpha_{1}$.

The numerical simulations indicate furthermore that, in the presence of the auto-amplified event $\mathcal{R}_{2}$, the division rates $\alpha_{2}$ and $\gamma$ have very little influence on $\mathbb{E}\left[\tau\right]$ while they affect dramatically $\mathbb{E}\left[C_{\tau}\right]$. It is also clear from the analytical solutions of the linear model, that the initial number of precursor cells $f_{0}$ and the spontaneous transition rate $\alpha_{1}$ have a major impact on both $\mathbb{E}\left[\tau\right]$ and $\mathbb{E}\left[C_{\tau}\right]$. A high sensitivity of model outputs to parameters is interesting to suggest possible key biological measurements (if feasible) in order to improve parameter identifiability (see paragraph 4.4.2-4.4.3).

4 Parameter calibration

In this section, we calibrate the model parameters using a likelihood approach. We first describe in subsection 4.1 the available experimental dataset, as well as in silico datasets that we use as a benchmark for our methodology. In subsection 4.2 we derive a likelihood function based on the embedded Markov chain from the underlying continuous-time Markov process. We explain how this likelihood is specifically adapted to the data, which are time-free measurements of cell numbers. We present the estimation results in subsection 4.3 for model ($\mathcal{M}_{FC}$) and each of the five different submodels: $(\mathcal{R}_{1},\mathcal{R}_{3})$, $(\mathcal{R}_{1},\mathcal{R}_{4})$, $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3})$, $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{4})$ and $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$. We recall that the different submodels are named by the reactions which have corresponding positive reaction rates. All the submodels considered are thus nested models, or reduced model compared to the full model $\eqref{Model_FC}$. We carry out a comprehensive comparison between the different models using model selection criteria. Thanks to a practical parameter identifiability analysis, we obtain model predictions in subsection 4.4, where we manage to retrieve hidden kinetic information and assess transit times and number of division events during the activation phase, with given confidence intervals. Finally, we discuss the biological interpretation of the calibration results in subsection 4.5.

4.1 Dataset description

Experimental dataset

Follicles undergoing the activation process have been classified according to three types braw-tal_studies_1997 ; gougeon_morphometric_1987 ; lundy_populations_1999 ; meredith_classification_2000 . Primordial follicles (Type I or B) have either not yet or just initiated activation; they are composed of a single layer of flattened cells surrounding the oocyte. Primary follicles (Type II or C) have completed initiation; they only contain cuboidal (transitioned) somatic cells organized in less than two layers (this means that some follicles are strictly mono-layered, while in others an extra partially fulled layer is being built-up). In between Types I and II lies a class of transitory follicles (Type IA or B/C), with a mixture of flattened and cuboidal cells coexisting within a single layer. The progression from Type I to Type II is accompanied with a more or less pronounced increase in the total cell number (flattened plus cuboidal cells) and enlargement in the oocyte (and follicle) diameter (see bottom-right panel of Figure 4).

We have made use of a dataset acquired in sheep fetuses lundy_populations_1999 ; wilson_01 (courtesy of Ken McNatty), which provides us with flattened (precursor) and cuboidal (proliferative) cell numbers in a sample of follicles distributed into the three activation steps. The dataset is subdivided into two subsets corresponding to two different sheep strains : the “wild-type” Romney strain and the “mutant” Booroola strain. The latter is characterized by a natural mutation affecting the receptor to growth factor BMP15 and resulting in the alteration of follicle development (see the Introduction section).

We have $90$ data points for the Wild-Type dataset, and $81$ for the Mutant dataset. More specifically, the measures consist of the cell numbers counted on the largest 2D cross-section of histologically fixed follicles of type I, IB or II. This 2D number can be correlated with the total 3D cell number from standard stereological considerations lundy_populations_1999 .

On a time horizon of several weeks, as it is the case in the experimental study we based on lundy_populations_1999 , a primordial follicle can undergo three different fates: it can either get activated, die, or remain quiescent. When considering the long term evolution (during the whole reproductive lifespan), all healthy primordial follicles will eventually get activated , the vast majority of non healthy primordial follicles will die before leaving the pool, and the remaining ones will never get activated reddy_10 (in women there are approximately one thousand follicles left in the ovaries after menopause – on the order of $0.2\%$ of the initial pool broekmans_09 ).

Except in the case of morphological abnormalities, it is not obvious to classify primordial follicles as healthy, in the sense that they will get activated at one time or another. To use appropriate data for the model fitting step, which should concern only, or at least mostly, “activable” follicles, we thus needed to gather complementary a priori biological information.

The viability of the oocyte enclosed in a primordial follicle is the main determinant of the follicle health. Whether an oocyte is viable or not at this developmental stage results from the process of follicle formation, that we describe here very briefly (we refer the interested reader to juengel_02 ; monniaux_18b ; tingen_09 ; wang_17 for a complete overview). During embryonic development, primordial germ cells colonize the territory corresponding to the future gonads. In females, these cells undergo several rounds of mitotic divisions while they interact locally with somatic cells. They gather in syncytium structures, the germ cell cysts, that ultimately fragment into primordial follicles. Only a small proportion of the germ cells survive to this step (25% in sheep smith_93 ), which also coincides with the entry into meiosis. Most oocytes will die and transfer a part of their cytoplasmic material to surviving oocytes. Each surviving oocyte recruits a variable number of somatic cells to build a primordial follicle sawyer_02 . The recruitment of a sufficient number of somatic cells is crucial to ensure future oocyte survival. Somatic cells secrete trophic factors, such as neurotrophin, whose intrafollicular levels should be high enough to guarantee the oocyte survival spears_03 . A further survival requirement is the existence of tight intercellular communications between the oocyte and its surrounding somatic cells, and between the somatic cells.

From these information, the following “viability criteria” can be derived to select proper type I follicles. First, there is a threshold oocyte diameter compatible with oocyte viability, and, second, for a given oocyte diameter, there is a threshold number of somatic cells surrounding the oocyte, that should be organized as a connected cell network paving the oocyte surface. The critical oocyte diameter can be determined in a rather straightforward way by comparison with the minimal oocyte diameter observed in transitory follicles. In the sheep species, following the study of cahill_81 (see also picton_01 for other species), we set the oocyte diameter threshold to 24 $\mu$m (corresponding to a follicle diameter threshold of 34 $\mu$m given the nominal depth of the somatic cell layer, and an absolute number of 15 cells in the largest 2D cross-section).

To assess the critical cell number relative to oocyte diameter, we computed a paving index, $e_{O}$, that represents the average contact length between a somatic cell and the oocyte: $e_{O}=\frac{\pi d_{O}}{N^{s}_{g}}$, where $d_{O}$ is the oocyte diameter and $N^{s}_{g}$ is the number of cells counted on the largest section. With a 24 $\mu$m diameter and 15 cells, we get a higher bound of $5\mu$m for $e_{O}$, that we applied as a filter to the rough data.

On the other side of the activation process, we have only retained the strictly mono-layered type II follicles. Indeed, we intend to deal with a final cell number as close as possible to the number reached at the first time when all flattened cells have transitioned to cuboidal cells (hence to the extinction time in the model), Yet, due to the oocyte enlargement and the resulting increased capacity of the first layer, one cannot preclude that a significant amount of cuboidal cells have been generated after the end of the transition period.

Combining these criteria, we get the dataset described on Figure 4, which illustrates the repartition of the data points according to the follicle type and sheep strain in each phase plane ($C$, $F$) ($C$, $p_{C}$), ($N$, $p_{C}$).

In silico datasets

In addition to the experimental dataset, we have constructed in silico datasets from simulations of SDE (1) in a way that mimics the experimental protocol (see details in Appendix 6.5).

In the sequel these datasets will be used as benchmark tools for the parameter identifiability study and the statistical comparison between the submodels and full model. In any case, the set of estimated parameters will match the set of cell events included in the model used to generate the in silico dataset. For instance, we will estimate the values of parameters $\alpha_{2}$ and $\gamma$ on the two datasets generated from submodel $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$ ($\alpha_{1}$ will be fixed to $1$ in the sequel).

4.2 Likelihood method

Since the experimental dataset is made of time-free observations, we are going to confront the model to the data using only the information on some state space values taken by the process, without their corresponding time information. This notion is intrinsically related to the embedded Markov chain which we detail below. We will use this Markov chain to compute a likelihood function. Note that the proliferative cell population increases by one cell at each event ($\mathcal{R}_{1}$, $\mathcal{R}_{2}$, $\mathcal{R}_{3}$ or $\mathcal{R}_{4}$), while the precursor cell population can either remain constant ($\mathcal{R}_{3}$ or $\mathcal{R}_{4}$) or decrease by one ($\mathcal{R}_{1}$ or $\mathcal{R}_{2}$). The proliferative cell population $C$ can thus be used as an event counter. Indeed, as a continuous-time Markov process, $X$ (defined in Eq. 1) can be decomposed into an embedded Markov chain $(F_{n},C_{n})_{n\in\mathbb{N}}$ and a sequence of random jump times $(s_{n})_{n\in\mathbb{N}}$ with

[TABLE]

Note that the sequence of jump times $(s_{n})_{n\in\mathbb{N}}$ corresponds exactly to the sequence of jump times associated with process C, and

[TABLE]

Given that $C_{n}=n$ is deterministic, it is clear that the precursor cell population $F_{n}$ (alone) is also a (non-homogeneous) Markov chain. To clarify the link with the data, we will index the embedded chain $F_{n}$ by the number of proliferative cells $C=c$, rather than by the number of events that occurred: let $F_{c}$ be the random variable corresponding to the number of precursor cells given that there are $c\in\mathbb{N}$ proliferative cells. According to the dichotomy between the two division events ($\mathcal{R}_{3}$, $\mathcal{R}_{4}$) and the two transition events ($\mathcal{R}_{1}$, $\mathcal{R}_{2}$), we deduce the law of $F_{c}$ at the “pseudo-time” $C=c$ from the law of $F_{c-1}$ at the “pseudo-time” $C=c-1$ as follows: for all $(f,c)\in\mathcal{S}$,

[TABLE]

where

[TABLE]

Hence $(F_{c})_{c\in\mathbb{N}}$ is a non-homogeneous discrete time Markov chain. Notice that the law of $C_{\tau}$, the number of proliferative cells at the extinction time of the precursor cells, corresponds to the law of the first “pseudo-time” $c$ such that $F_{c}=0$, e.g. $C_{\tau}=\inf\{c\in\mathbb{N}^{*},\,F_{c}=0\}$.

In addition to Eq. (15), to compute the law of $(F_{c})$, we need to specify an initial condition $F_{0}$. As detailed in Section 4.1, the initial pool of flattened cells is highly variable. To limit the number of parameters, we assume that the initial number of precursor cells follows a truncated Poisson law on $\mathbb{N}^{*}$ (with a single parameter $\mu\in\mathbb{R}_{+}$) given by, for all $f\in\mathbb{N}^{*}$,

[TABLE]

Then, we can use Eq. (15) to compute $\mathbb{P}[F_{c}=f]$ by recurrence from the initial probability vector $(\mathbb{P}[F_{0}=i])_{i\in\llbracket 0,c+f\rrbracket}$. Hence, we have built a discrete time Markov chain $(F_{c})_{c\in\mathbb{N}}$ from model ($\mathcal{M}_{FC}$) adapted to our time-free observations.

As can be seen from Eq. (16), the timescale cannot be inferred, so that we fix arbitrarily $\alpha_{1}=1$, whatever the dataset, to obtain dimensionless parameters. The time unit of the remaining parameters is thus relative to the timescale of one spontaneous transition event, and their estimated values may depend on the specific dataset (experimental or in silico).

Finally, we suppose that all data points are independent of one another, and that the observations are free of measurement errors. Therefore, accordingly to Eq. (15)-(16)-(17), our statistical model assumes that the observed variability is due to a random initial number of precursor cells, and to the occurrences of random cell events among cell transitions and cell divisions. We obtain the following likelihood function

[TABLE]

for $N$ data points $(f_{i},c_{i})_{i=1..N}$ and the parameter vector $\theta$ depending on each submodel.

We infer the parameter values using the maximum likelihood estimator (MLE), and apply the practical approach based on profile likelihood estimate (PLE) to analyze the parameter identifiability and assess confidence intervals raue_structural_2009 . We also perform model selection using classical AIC and BIC criteria to discriminate between the full model and different submodels. The whole procedure is described in Appendix 6.6.

Note that the initial condition parameter $\mu$ can be either estimated together with the other parameters from a given dataset, using the likelihood given in Eq. (18), or, alternatively, from the cell number of the primordial follicles only. In the latter case, with the law of $F_{0}$ given by (17), we obtain the likelihood function

[TABLE]

for $N^{\prime}$ data points $(f_{i},0)_{i=1..N^{\prime}}$. From the likelihood defined in Eq. (19), we deduce MLE and PLE to infer the value of $\mu$ solely from the primordial follicle data.

4.3 Fitting results

In this subsection, we present our fitting results using the procedure described in subsection 4.2 for several submodels derived from model ($\mathcal{M}_{FC}$):

• •

two-event submodels, including the spontaneous transition event together with either the asymmetric $(\mathcal{R}_{1},\mathcal{R}_{3})$ or symmetric division $(\mathcal{R}_{1},\mathcal{R}_{4})$;

• •

three-event submodels, the linear submodel $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$ and the two-nonlinear one, including auto-amplified transition events, together with either the asymmetric $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3})$ or symmetric $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{4})$ division event;

• •

the full model ($\mathcal{M}_{FC}$)$=(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3},\mathcal{R}_{4})$

The fitting results obtained with the total likelihood (Eqs. (15)-18) on the experimental datasets are shown in Figure 5 for submodels $(\mathcal{R}_{1},\mathcal{R}_{3})$, $(\mathcal{R}_{1},\mathcal{R}_{4})$ and the full model . The corresponding fitting results for the in silico datasets for the same submodels are provided in Figure 11 (Appendix 6.7). The fitting results for the three-event submodels for both the experimental and in silico datasets are provided in Appendix 6.7. For both the Wild-Type and Mutant datasets, a visual inspection shows that submodel $(\mathcal{R}_{1},\mathcal{R}_{4})$ leads to a “direct” transition, followed by prolonged cell proliferation after precursor cell extinction, while, with submodel $(\mathcal{R}_{1},\mathcal{R}_{3})$, there is a higher probability that the total number of cells increases before precursor cell extinction. The model selection criteria, summarized in Table 1, shows that all submodels without cell event $\mathcal{R}_{4}$ can be safely rejected. The visual inspection of Figure 5 leads to the following explanation. If event $\mathcal{R}_{4}$ is present, as in submodel $(\mathcal{R}_{1},\mathcal{R}_{4})$, the proliferative cells can keep dividing after the extinction of the precursor cells (line $F=0$). Once the precursor cell number reaches zero for a given $c$, all remaining points $(0,c^{\prime})$ for $c^{\prime}\geq c$ are reached with probability one, which results in a high contribution of all $(0,c)$ data points to the maximum likelihood. In contrast, if event $\mathcal{R}_{4}$ is not present, as in submodel $(\mathcal{R}_{1},\mathcal{R}_{3})$, the process stops as soon as the precursor cell population $F$ gets extinct, which prevents the likelihood of all $(0,c^{\prime})$ points from being close to one (they rather take all intermediate values). This observation is consistent with the fitting results of the in silico datasets (Figure 10 in Appendix 6.7).

The model selection criteria further suggest that the best models associated with the experimental datasets are the full model and the three-event linear submodel $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$. The two-event submodel $(\mathcal{R}_{1},\mathcal{R}_{4})$ appears to be a possible alternative but still less relevant than the two others. In Figure 5, we observe that the trajectories associated with an intermediate level of cell proliferation before precursor cell extinction are more likely in the full model than in the two-event submodel $(\mathcal{R}_{1},\mathcal{R}_{4})$, with a more pronounced effect for the mutant subset than the wild-type subset. We will come back to this consideration in section 4.4.

The parameter identifiability study (detailed in Appendix 6.7) leads to the following results:

• •

The initial condition parameter $\mu$ is always practically identifiable (Figures 6 and 12) and the MLE yields similar values from one submodel to another (see Tables 4 and 5). Moreover, its fitted value is close to the true one for the in silico datasets, with some small bias in some cases (Figures 11 and 12).

• •

The asymmetric division parameter $\alpha_{2}$ is tightly identifiable only in submodel $(\mathcal{R}_{1},\mathcal{R}_{3})$, (Figure 6). However, significant bias occurs, as revealed by the in silico datasets yet (Figures 11 and 12).

• •

The symmetric division rate $\gamma$ is rarely practically identifiable, yet an upper-bound can always be obtained (Figures 6, 11 and 12).

• •

The self-amplified transition rate $\beta$ is not identifiable in any case (Figures 6, 11 and 12).

This one-dimensional parameter identifiability analysis hides however more subtle parameter constraints. The self-amplification transition rate is actually constrained to be greater than the symmetric division rate $\gamma$, as shown in the two-dimensional profile likelihood analysis in Figure 13 in Appendix 6.7. This result confirms the tendency observed with the best fit trajectories in Figure 5, that favor transition over proliferation.

The fitting results obtained on the initial condition parameter $\mu$ from primordial follicle data (using the likelihood given by Eq. (19)) is shown in Figure 7. We have obtained identifiable parameter values with each submodel, yet associated with broader confidence intervals than with the global fitting approach given by Eqs.(15)-(18). As expected, using more information reduces the uncertainty, hence the confidence intervals are smaller when the whole datasets are used (for all models and subsets considered).

4.4 Model prediction

In this subsection, we use the MLE together with their confidence interval obtained with the PLE of the best models (the two linear submodels $(\mathcal{R}_{1},\mathcal{R}_{4})$ and $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$ and the full model) to infer information on the experimental subsets.

4.4.1 Distribution of the initial condition

In Figure 7, we compare for both the Wild-Type and Mutant subset the distributions derived from model $(\mathcal{R}_{1},\mathcal{R}_{4})$, $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$ and $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3},\mathcal{R}_{4})$, using the whole data, together with the distribution inferred from the primordial follicle data only. From the top panels of Figure 7, we observe that in all cases, there is an overestimation of the head and tail of the distribution of $F_{0}$, which suggests that a more peaked distribution than the truncated Poisson distribution would be more suitable for the initial condition. The distribution inferred from the primordial follicle data only is slightly closer to the datapoint than the distribution with $\mu$ inferred using the complete follicle data (as expected), as assessed by the evaluation of the likelihood (19) at each MLE, shown in the lower panels of Figure 7.

A detailed inspection of the lower panels of Figure 7 shows furthermore that the likelihood (19) based on the primordial follicle data cannot discriminate between the Wild-Type and Mutant subset. However, using the likelihood (15)-(18) with the whole data induces a shift of approximately one cell in average, in opposite directions for the Wild-Type and Mutant subset: for the Wild-Type subset, the mean cell number is found to be greater when the whole data are used, while for the Mutant subset, the mean cell number is found to be smaller (for all three models considered). Hence, considering the subsequent follicle trajectories, shaped by transition and proliferation, modifies the most likely value of $\mu$ and can discriminate the Wild-Type subset from the Mutant subset. The precise value of $\mu$ is biologically important, since it can be considered as the equivalent of the number of founder cells in lineage studies. Indeed, until ovulation (where the total cell number is on the order of several millions in sheep), there will not be any recruitment of somatic cells, and all cells with derive from the initial flattened cells.

4.4.2 Proliferative cell proportion: reconstruction of time

In Figure 8, we represent the predicted changes in the proliferative cell proportion with respect to time. For sake of readability, these predictions are derived from the deterministic formulation of the full model (Eq. (22)). We expect that a similar trend would be observed with the stochastic CTMC formulation. For each model, we superimpose the time trajectories corresponding to the parameter combinations for which the PLE is below the $95\%$ threshold. In both the Wild-Type and Mutant cases, despite the uncertainty affecting the model parameters for the two linear submodels (left and right upper panels), the dynamics just exhibit small uncertainties: the proportion of proliferative cells reaches $50\%$-$70\%$ in one time unit, which corresponds to the time unit of a single spontaneous transition event. This might due partly to the fact that parameter $\gamma$ is partially identifiable and is estimated to relatively low values. In contrast, the lack of parameter identifiability of the full model results in a huge uncertainty on the dynamics, that can be up to 5 order of magnitude faster than a single spontaneous transition event: the proportion of proliferative cells reaches $50\%$ between $10^{-6}$ and $1$ time unit. Indeed, cell event $\mathcal{R}_{2}$ (controlled by parameter $\beta$) can speed up the transition dynamics, and cell event $\mathcal{R}_{3}$ (controlled by parameter $\alpha_{2}$) can trigger the first transition, leading to a possible fast activation which avoids the bottleneck of the spontaneous transition timescale ($\alpha_{1}=1$). It is difficult to instantiate these relative durations in physical time units. The only kinetic information available on the activation process is given by studies that have monitored the sequential apparition of different follicle types during fetal development. In wild-type animals, the first primordial follicles appear around 75 days of gestation, while the first primary follicles are observed around 100 days mcnatty_development_1995 . A 25 day-duration can thus be considered as close to the minimal duration. No clear timescale separation between the Wild-type and Mutant dynamics can be revealed, although some parameter combinations are compatible with a faster transition in the Wild-Type case than in the Mutant case. This is again compatible with monitoring studies, which observed that the times of apparition of both the first primordial and primary follicles are shifted compared with wild-type animals (they appear a little later), yet the delay in between does not appear to be significantly different.

4.4.3 Mean extinction time, mean number of cells at the extinction time and mean number of division events during activation

In Figure 9, we represent the mean number of proliferative cells, $\mathbb{E}\left[C_{\tau}\right]$, and the mean number of division events during activation, $\mathbb{E}\left[C_{\tau}-F_{0}\right]$, as a function of the mean extinction time $\mathbb{E}\left[\tau\right]$, as predicted from the selected (sub)models $(\mathcal{R}_{1},\mathcal{R}_{4})$, $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$ and $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3},\mathcal{R}_{4})$. These predictions are obtained from a direct stochastic simulation of the trajectories of each model (with Gillespie algorithm, or SSA)333We use here the direct simulation rather than Algorithm 1, because the parameter range explored by the symmetric division rate $\gamma$ gets close to the theoretical necessary and sufficient condition $\gamma<\alpha_{1}+\beta$, while the Algorithm 1 requires $2\gamma<\alpha_{1}+\beta$., using the parameter values obtained from the identifiability analysis, for which the PLE is below the $95\%$ threshold. For each subset (Wild-Type or Mutant), the predicted value for $\mathbb{E}\left[C_{\tau}\right]$ is similar in each submodels and lies between 8 and 10 cells. Interestingly, the predicted value for $\mathbb{E}\left[C_{\tau}\right]$ is approximately 6-8 cells lower than the empirical mean number of proliferative cells obtained directly from the primary follicle data (data points $(0,C)$ with $F=0$) (Figure 9, top panels). This observation is consistent with the trajectory analysis performed from Figure 5, from which we have concluded that the activation process follows with high probability a trajectory reaching state $F=0$ with a low cell number, and characterized by direct transition and very little concomitant cell proliferation. Similarly, $\mathbb{E}\left[C_{\tau}-F_{0}\right]$ is approximately 5-7 cells lower than the increase in the mean empirical number of cells between the primordial follicle datasets and primary follicle datasets (Figure 9, bottom panels). $\mathbb{E}\left[\tau\right]$ in the two linear submodels $(\mathcal{R}_{1},\mathcal{R}_{4})$ and $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$ depends only on the initial condition and is estimated to a value around $2.5$ a.u. with a small uncertainty, similarly as in Figure 8. In contrast, the full model yields a larger uncertainty on $\mathbb{E}\left[\tau\right]$, with a confidence interval between $10^{-6}$ and $0.5$ a.u. for the Wild-Type subset, and between $10^{-6}$ and $2.5$ a.u. for the Mutant subset, consistently with the prediction on the dynamics of the proliferative cell proportion (Figure 8). From our theoretical results on parameter sensitivity in section 3.3 (see Figure 3), we have found that $\beta$ has a profound impact on $\tau$. Any additional knowledge on the follicle activation duration would thus be valuable to further constraint the parameter uncertainty.

Predictions on the mean number of divisions events could be in theory amenable to validation by experimental cell kinetics study. Such studies, enabling for instance to infer the possibly time-varying doubling times in cell populations have been performed for later developmental stages (turnbull_77 in sheep or pedersen_70 in mice). They cannot be conducted as such for the earliest stages because of the excessive slowness of cell events. Promising ex-vivo/in vitro devices morohaku_16 reproducing all steps of follicle development could be appropriate to settle elaborate cell lineage tracing informing on cell division events. Yet this is a long term perspective, since such devices are rather at the proof-of-concept level for the time being. In addition, they are up to now restricted to the mouse, since no feasible culture system for primordial and primary follicles is yet available for other species (see the overview picture including sheep in morohaku_19 ).

4.5 Biological interpretation

From the primordial follicle data, we have found that the mean initial number of precursor cells for the Wild-Type subset is about the same as for the Mutant. Moreover, the prediction on the total number of proliferative cells at the end of the activation phase, $\mathbb{E}\left[C_{\tau}\right]$, is also very similar in the Wild-Type and Mutant cases. The observed shift in opposite directions for the mean initial cell number inferred from the MLE of the dynamical models (see bottom panel of Figure 7) is thus compensated for by the differences in cell dynamics. The number of divisions during the transition is smaller in the Wild-Type than in the Mutant subset ($\mathbb{E}[C_{\tau}-F_{0}]\approx 2$ in Wild-Type, $\mathbb{E}[C_{\tau}-F_{0}]\approx 4$ in Mutant), as a result of a global difference between the MLE parameters: the order of magnitude of the division rates are closer to that of the transition rates in the Mutant compared to the Wild-Type subset. In overall, we conclude from our extensive datafitting analysis that the Wild-Type subset exhibits a clearer separation of dynamics during follicle activation (first cell transition, then cell proliferation), while in the Mutant cell proliferation could occur at a substantial rate before precursor cell extinction. We note that this conclusion has to be tempered by the sparse character of our experimental dataset. In particular, a detailed examination of the experimental data reveals that the four data points available for transitory follicles in the Wild-Type subset correspond to a clearly higher number of precursor cells than any of the primordial follicles, which certainly impacts our results. In contrast, the Mutant subset contains transitory follicles with significantly fewer precursor cells than the primary follicles.

Even if there is a clear trend in the data to substantiate the existence of an auto-amplification of the transition from flattened to cuboidal cells, complementary data would be useful to decide the question. Indeed, getting more datapoints with a proportion of cuboidal cells in the range of 50 to 100 % would constrain a step further the follicle activation trajectories, hence the parameter values and differences between nested models. The very fact that fewer follicles are counted in this range in the Wild-Type subset pleads for a possible acceleration of the transition. Statistically, including more animals and more gestation times in the study would increase the number of data, including data missing in our current dataset, yet it would require enrolling many experimental animals. Ideally, monitoring the cell dynamics of ovarian follicles in vivo, in a non invasive manner, would provide all needed data. Yet, it is far from being a reachable target at the moment. Even the morphological monitoring of follicles (individual changes in follicle diameters) can only be performed for much later developmental stages due to size resolution (no reliable data can be obtained below 2mm diameter). An alternative would be to record the location of the cuboidal cells with respect to the flattened ones, in consistency with the spatial interpretation of auto-amplification. The auto-amplification rate is motivated by two possible (and non exclusive) underlying mechanisms. First, the very first cell transitions could awake the oocyte and settle a positive feedback loop between the somatic cells and the oocyte adhikari_09 ; monniaux_18b that would in turn secrete stimulatory factors reaching the surrounding somatic cells by diffusion (global amplification). Second, communications between adjacent somatic cells could help propagate activation step by step, from one (or a few) originally activated cell (local amplification). Local amplification might be detected in the data by recording the location of cuboidal cells and checking whether cluster of spatially related cuboidal cells can be detected. Global detection is expected to have a more homogenous effect, hence to be hardly detectable from static histological data.

Finally, we highlight that the $\beta$-free linear submodel $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$ performs as well as, and even better than the complete model ($\mathcal{M}_{FC}$) $(\mathcal{R}_{i})_{i\in\llbracket 1,4\rrbracket}$ in Mutant compared to Wild-Type ewes, which is compatible with the functional hypotheses applicable to the BMP15R mutation reader_12 . Indeed, one could speculate that the diminished BMP15 signaling would hamper the molecular dialog between the oocyte and somatic cells after follicle activation triggering, so that the auto-amplified cell event would barely occur in the Mutant group.

5 Conclusion

In this work, we have introduced a stochastic nonlinear cell population model to study the sequence of events occurring just after the initiation of follicle growth. We have characterized the dynamics of precursor and proliferative cell populations according to the parameter values, for both the stochastic model and its deterministic mean-field counterpart. We have studied in details the extinction time of the precursor cell population, and designed an algorithm to compute numerically both the mean extinction time and mean number of proliferative cells at the extinction time. The algorithm is based on a domain truncation similar to the Finite State Projection (FSP) method proposed in munsky_finite_2006 ; kuntz_deterministic_2017 . The FSP approach aims to approximate the law of the process at a given time by solving a truncated version of the Kolmogorov forward system. We have adapted the FSP algorithm to close the infinite recurrence relation satisfied by the extinction time moments. We have found a consistent spatial boundary to solve the closure problem, thanks to a coupling technique and tractable upper-bound process. The numerical cost of the algorithm is deeply related to the proper choice of the upper-bound processes and gets worse than direct simulation as $2\gamma$ gets close to the required bound $\alpha_{1}+\beta$ of Algorithm 1.

This algorithm has nevertheless allowed us to investigate the parameter influence on the precursor cell extinction time and number of proliferative cells at the end of the follicle activation phase. The auto-amplified transition rate $\beta$ exerts a critical control on the mean extinction time, with a sharp timescale reduction when $\beta$ exceeds the spontaneous cell transition $\alpha_{1}$, while the division rates ($\alpha_{2}$, $\gamma$) have relatively less effect. The effect of the auto-amplification process is probably dependent on the specific parameterization of the cell event rates chosen in this work, yet our findings bring interesting insight into the mechanisms underlying follicle activation; nonlinear feedbacks mediated through cell-to-cell communication certainly play a role, and our estimation results have shown that any impairment of this feedback would change drastically the kinetics of follicle activation.

Moreover, our results can be useful to understand the variability in the cell numbers among ovarian follicles at the end of the activation phase, which can be used as initial conditions for models describing the following stages of follicle development clement_coupled_2013 ; CRY2019 . Going even further, the sequence of events occurring just after the initiation of follicle growth is determinant for the remaining of the entire follicle development process. The whole cell population in mature (ovulatory) follicles (up to tens of millions in large mammal species as humans) emanates from the few cells a primordial follicle is endowed with. The timings of the cells’ first divisions will determine the distribution of cytological cell ages in the population, which will ultimately influence the distribution of the times of cell cycle exit in fully differentiated cells. Collectively, the exit time distribution controls the switch from proliferation to differentiation, a key event in the selection of ovulatory follicles clement_13b . Also, the proliferative vitality of the cuboidal (transitioned) cells will control the clonal composition of the follicles and participate in the spatial and functional heterogeneity within follicle cell populations, persisting very late in development.

We have performed the parameter calibration in a special context of time-free data. It turns out that the proliferative cell number can be seen as a clock for the whole process, and that the embedded Markov chain is better adapted to time-free data than the continuous-time model. We have used the embedded Markov chain to define a proper likelihood function and a statistically rigorous framework. The likelihood function has allowed us to perform an extensive data fitting analysis, using the very useful concept of profile likelihood estimate. This analysis sheds light onto several aspects of the activation of ovarian follicles. First, the transition scenario, where cell proliferation is mostly posterior to cell transition, and the cell number increase is moderate, seems to be predominant versus a more proliferative scenario. While the question is still open, it seems likely that cell transition is favored in the Wild-Type strain compared to the Booroola mutant strain. With the available experimental dataset, we have yet not managed to make a clear distinction between, on one side, a progressive transition with a steady net flux from flattened to cuboidal cells, and, on the other side, an auto-catalytic transition with an ever increasing flux all along the activation phase.

Beyond our application in female reproductive biology, we believe that the modeling approach presented here can have a more generic interest in cell kinetics related issues, especially when a small number of cells is involved. Also, from the mathematical biology viewpoint, the analysis performed on the extinction time, combining theoretical (coupling) and numerical (finite state projection) tools may have an interest for first passage time studies in stochastic processes.

6 Appendix

6.1 Justification of the choice of the rate of $\mathcal{R}_{2}$

As detailed in Section 4.5 the auto-amplification can result from two non-exclusive mechanisms, a nonlocal (global) one and a local one.

Global amplification: consider that each proliferative cell sends a fixed amount of growth signals to the oocyte. The oocyte thus receives a signal proportional to the number of proliferative cells $C$. We consider that the oocyte secrete in turn (instantaneously) a stimulatory signal, at a level proportional to the amount of growth signals received from somatic cells. By homogeneous diffusion, the oocyte signal is shared equally to all somatics cells, so that each precursor cell receive a signal proportional to $C/(F+C)$.

Local amplification: for a given precursor cell, assuming a random repartition of the cell types around the oocyte (hence neglecting local cell-to-cell effects), the probability that a neighbor cell is a proliferative cells is $C/(F+C-1)$, which is also consistent with our choice.

6.2 Mean-field formulation

To get some insight into the model behavior, we describe the mean-field version of model $\mathcal{M}_{FC}$, given by the following set of ODE:

[TABLE]

with the initial condition $(f(0),c(0))=(f_{0},0)$, with $f_{0}\in\mathbb{R}_{+}$. We start by solving analytically the deterministic formulation, and then investigate the effect of each parameter on the model outputs.

From the ODE sytem (20), we deduce the change in the proliferative cell proportion $p_{C}(t):=\frac{c(t)}{f(t)+c(t)}$:

[TABLE]

From ODEs (20) and (21), using the classical method of separation of variables, we can compute the analytical expressions for the proliferative cell proportion $p_{C}(t)$, proliferative cell number $c(t)$ and precursor cell number $f(t)$:

Proposition 5

The solution of the ODE system (20) is, for all $t\geq 0$,

[TABLE]

In addition, the solution of ODE (21) is

[TABLE]

and the total cell number verifies

[TABLE]

From Proposition 5, it is clear that the proliferative cell proportion $p_{C}$ converges to $1$. If $\gamma>0$, the proliferative cell number $c$ grows asymptotically exponentially at a rate $\gamma$ when $t\to\infty$. If $\gamma=0$, $c(t)$ is bounded because $t\mapsto 1-p_{C}(t)$ is converging exponentially fast to [math], hence is integrable on $(0,\infty)$. Moreover, the proliferative cell proportion $p_{C}$ has an inflexion point if and only if

[TABLE]

An inflexion point denotes the presence of at least two distinct phases, with a first progressive acceleration phase followed by a saturating phase.

Finally, note that according to the observed variables, the submodels cannot be distinguished from one another, or, alternatively, different parameter values (within a same submodel) may lead to identical outputs. Indeed, the changes in the precursor cell population are independent of parameters $\alpha_{2},\gamma$, and, more strikingly, parameters $\beta$ and $\gamma$ cannot be separated in the analytical solution (22), leading to the same kinetic patterns for $p_{C}$ as long as the combination $\gamma+\beta$ remains unchanged.

6.3 Analytical expressions in the linear case

Proof (Proof of Proposition 1)

Let $t\geq 0$ and $f\in\llbracket 0,f_{0}\rrbracket$. Since $F_{t}$ is autonomous and is a pure death process, we can directly write the following forward Kolmogorov equation: for all $f\in\llbracket 0,f_{0}\rrbracket$,

[TABLE]

Solving by recurrence (23), we deduce that, for all $f\in\llbracket 0,f_{0}\rrbracket$,

[TABLE]

Note that $\mathbb{P}\left[F_{t}^{L}=0|F_{0}=f_{0}\right]=(1-e^{-\alpha_{1}t})^{f_{0}}$ which converges to $1$ when $t\mapsto 1-p_{C}(t)$. Hence, process $F^{L}$ extincts almost surely (a.s.) when $t$ goes to infinity, hence $\tau_{L}<\infty$. Before computing the law of $\tau_{L}$, we can directly obtain its mean using the recursive expression (4):

[TABLE]

Using again Eq. (4), we deduce that $\tau_{L}(=T_{f_{0}})$ follows a generalized Erlang law whose density function is:

[TABLE]

Due to the specific form of the exponential rate, we can simplify Eq. (24) further. As $\displaystyle\prod_{j\neq i,j=0}^{f_{0}-1}(f_{0}-j)=\frac{f_{0}!}{f_{0}-i}$ and

[TABLE]

we deduce

[TABLE]

$\square$

Proof (Proof of Proposition 2)

According to Proposition 1, $\tau_{L}$ is a.s. finite. To take the expectation of $C^{L}_{t}$ at time $t=\tau_{L}$, we check that $\mathbb{E}\left[C^{k,j}_{\tau_{L}-T_{k}^{j}}\right]<\infty$, for all $k$ and $j$. For all $t\geq 0$, $C_{t}^{k,j}$ is $L_{1}-$integrable (as a Yule process) with $\mathbb{E}\left[C^{k,j}_{t}\right]=e^{\gamma t}$. Conditionning on the law of $\tau_{L}$, we get (with the change of variables $x=1-e^{-\alpha_{1}t}$)

[TABLE]

where $B$ is the standard Beta function. Hence $I<\infty$ if and only if Hypothesis 3 holds. Note that using the properties of the Beta function, we have

[TABLE]

where we use the notation $\left(m-x\right)!=\prod_{k=1}^{m}(k-x)$.

Thus, if Hypothesis 3 holds true, and given that $C^{k,j}$ is a positive increasing process, we deduce:

[TABLE]

Then, taking the expectation of (6) at time $t=\tau_{L}$, we obtain:

[TABLE]

Moreover, we have that each counting process $N_{k}(t)$ can be dominated by

[TABLE]

so that

[TABLE]

Finally, conditionally on $\tau_{L}$, $\mathcal{Y}_{3}(\tau_{L})$ is independent of each $C^{k,j}_{\tau_{L}}$, and the latter are independent and identically distributed random variables. Using that

[TABLE]

and the Wald equation (Feller, , Chap. XII) , we obtain

[TABLE]

which is finite under Hypothesis 3. Finally, if Hypothesis 3 does not hold, we have, as long as $f_{0}\geq 2$:

[TABLE]

In some special cases, Formula (26) can be used to obtain the first moment of $C^{L}_{\tau_{L}}$.

When $\gamma$ is zero, then for all $t\geq 0$, for all $k\in\llbracket 1,f_{0}\rrbracket$ and for all $j\in\llbracket 1,N_{k}(\tau_{L})\rrbracket$, $C^{k,j}_{t}=1$. We deduce directly from Eq. (26) that

[TABLE]

From Eq. (7), we have

[TABLE]

by Poisson process property. Since for all $t\in[T_{k},T_{k+1})$, $F^{L}_{t}=f_{0}-k$, we deduce that $\mathbb{E}\left[N_{k}(\tau_{L})\right]=\mathbb{E}\left[\alpha_{2}(f_{0}-k)(T_{k+1}-T_{k})\right]$. Using (4), we deduce that $\mathbb{E}\left[N_{k}(\tau_{L})\right]=\frac{\alpha_{2}(f_{0}-k)}{\alpha_{1}(f_{0}-k)}=\frac{\alpha_{2}}{\alpha_{1}}$ and conclude with (27).

When $\alpha_{2}$ is zero, $N_{k}(t)$ is null for all $t\geq 0$, and we deduce directly from (26) that

[TABLE]

Since $T_{f_{0}}=\tau_{L}$, we have $C^{f_{0},0}_{\tau_{L}-T_{f_{0}}}=1$. Let $k\in\llbracket 1,f_{0}-1\rrbracket$. Since $\tau_{L}-T_{k}\overset{(law)}{=}\sum_{i=k+1}^{f_{0}}\mathcal{E}\left(\alpha_{1}(f_{0}-i+1)\right)\overset{(law)}{=}\sum_{i=1}^{f_{0}-k}\mathcal{E}\left(\alpha_{1}i\right)$, using Proposition 1, we deduce that the density function of $\tau_{L}-T_{k}$ is

[TABLE]

Then, conditioning $C^{k,0}_{\tau_{L}-T_{k}}$ on the law of $\tau_{L}-T_{k}$, we first deduce that

[TABLE]

Then, since $\mathbb{E}\left[C^{k,0}_{t}\right]=e^{\gamma t}$, we have, similarly as in Eq. (25),

[TABLE]

which ends the proof using (28).

$\square$

The following proposition is analogous to Proposition 2, yet with the decoupled processes $\tilde{F}$ and $\tilde{C}$, whose moments are easier to estimate. Note that parameters $\tilde{\alpha},\tilde{\beta},\tilde{\gamma}$ below are generic ones.

Proposition 6

Let $\tilde{F},\tilde{C}$ be independent pure-jump stochastic processes on $\mathbb{N}$, of infinitesimal generators

[TABLE]

with deterministic initial condition $\tilde{F}(0)=f_{0}$ and $\tilde{C}(0)=n\geq 1$, and where $\tilde{\alpha},\tilde{\beta},\tilde{\gamma}$ are non-negative rate parameters. Let

[TABLE]

For any $p\geq 1$,

[TABLE]

if, and only if,

[TABLE]

Moreover, we have:

• •

if $\tilde{\gamma}>0$: for $p=1$,

[TABLE]

and for $p=2$,

[TABLE]

• •

if $\tilde{\gamma}=0$:

[TABLE]

[TABLE]

Proof

Since $\tilde{\tau}$ and $\tilde{C}$ are independent, we deduce by conditioning on $\tilde{\tau}$ that

[TABLE]

where $f_{\tilde{\tau}}$ is the density probability of $\tilde{\tau}$. Since $\tilde{F}$ is linear, we apply Proposition 1 and obtain

[TABLE]

Now, we suppose that $\tilde{\gamma}>0$. Then, $\tilde{C}$ can be decomposed as the independent sum of $n$ Yule processes starting from $1$ (see Eq. (5)) and a birth process with immigration (starting from [math]). It is classical that the Yule process follows a geometric law of parameter $e^{-\tilde{\gamma}t}$, and the birth process with immigration follows a negative binomial law $\mathcal{BN}\left(\frac{\tilde{\beta}}{\tilde{\gamma}},e^{-\tilde{\gamma}t}\right)$, there exists $k,K>0$ (depending on model parameters, but independent of $t$) such that, for all $t\geq 0$,

[TABLE]

Combining Eq. (33) with Eqs. (31) and (32) yields (29). To obtain the remaining analytical formulas, we note that

[TABLE]

and

[TABLE]

Also, for any $p$ such that (30) holds true, we have (with the change of variables $x=1-e^{-\tilde{\alpha}t}$)

[TABLE]

where $B$ is the standard Beta function. We deduce that

[TABLE]

where we use the notation $\left(m-x\right)!=\prod_{k=1}^{m}(k-x)$. Then, using Eqs. (34)-(35) and (32), we deduce from (31) that

[TABLE]

and

[TABLE]

If $\tilde{\gamma}=0$, then $\tilde{C}$ is a pure immigration process starting from $n$, and follows a shifted Poisson law $n+\mathcal{P}\left(\tilde{\beta}t\right)$ at time $t\geq 0$. Using the same approach, we obtain that

[TABLE]

and

[TABLE]

.

$\square$

6.4 Numerical scheme for $\mathbb{E}[\tau]$ and $\mathbb{E}[C_{\tau}]$

Pseudo-code

We design algorithm 1 to compute a numerical estimate of $g(f_{0},0)$, solution of Eq. (13) that represents either $\mathbb{E}\left[\tau\right]$ or $\mathbb{E}\left[C_{\tau}\right]$ according to the specific choice of boundary condition. This algorithm requires $\gamma<\alpha_{1}+\beta$ to compute $\mathbb{E}\left[\tau\right]$, and $2\gamma<\alpha_{1}+\beta$ to compute $\mathbb{E}\left[C_{\tau}\right]$, in agreement with Theorem 3.1, Proposition 3 and Proposition 4. The prefactor $A$ given below is obtained thanks to Proposition 6.

6.5 In silico dataset

We generate in silico datasets to further explore parameter identifiability. For each submodel, we choose two different parameter sets with contrasted values in the division rates $\alpha_{2}$ or $\gamma$ and/or transition rate $\beta$. The parameter values are summarized in Table 2. We obtain the corresponding 10 datasets by simulating $1,000$ trajectories from the SDE (1), with the Gillespie algorithm gillespie_general_1976 , starting from the initial condition $(F_{0},0)$ at time $t=0$ up to the time when $C(t)=31$ (the value $C(t)=31$ corresponds to the maximal number of cuboidal cells observed in the experimental dataset). The initial random variable $F_{0}$ follows a truncated Poisson law of parameter $\mu$ (see Eq.(17)). For each trajectory, we select uniformly randomly one point $(f,c)$ among the state space points reached by the trajectory, so that each in silico datasets is composed of $N=1,000$ points. This way of sampling, letting to time-free and uncoupled datapoints, mimics the experimental protocol.

6.6 Detailed fitting procedure

Maximum likelihood estimator

For each submodel and dataset, the optimal parameter values are given by the MLE $\hat{\theta}=\left(\widehat{\beta},\widehat{\alpha_{2}},\widehat{\gamma},\widehat{\mu}\right)$, which we compute by minimizing the negative log-likelihood,

[TABLE]

for a dataset $\mathbf{x}$ and where $\Theta$ is constructed by fixing all parameters related to the nonpresent events to the singleton $\{0\}$: for instance, in submodel $(\mathcal{R}_{1},\mathcal{R}_{4})$, we have $\Theta=\{0\}\times\{0\}\times\mathbb{R}_{+}\times[1,+\infty)$.

To compute the minimum, we use a derivative-free optimization algorithm: the Differential Evolution (DE) algorithm storn_differential_1997 . In the following, we describe the whole procedure for the complete model $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3},\mathcal{R}_{4})$. The algorithm starts from an initial population in which each individual is represented by a set of real numbers $(\beta,\alpha_{2},\gamma,\mu)$. Then, the population evolves along successive generations by mutation and recombination processes. At each generation, the likelihood function is used to assess the fitness of the individuals, and only the best individuals are kept in the population. We have set the intrinsic optimization parameters as follows: the initial population has a size of 20 individuals, and the probability of mutation and crossing-over equals to 0.8 and 0.7 respectively. The starting individual parameter sets are defined on a log scale, and drawn from a uniform distribution on $\Theta=[-6,6]^{3}\times[0,1.5]$. The algorithm was run over 1,000 iterations.

Profile likelihood estimate

For each $i$th component of the MLE $\hat{\theta}_{i}$, $i\in\llbracket 1,4\rrbracket$, we compute a vector $\hat{\theta}|[\theta_{i}=x]$ on a grid $G_{i}$ around the MLE $\hat{\theta}$, with $x\in G_{i}$:

[TABLE]

and its associated PLE (vector) $\mathcal{L}(\mathbf{x};\hat{\theta}|\theta_{i})$. We design the grid $G_{i}$ around the MLE $\hat{\theta}_{i}$ with a fixed step size (see Table 3 for details), and re-optimize the remaining parameters using the DE algorithm with the same optimization parameters (mut=0.8, crossp=0.7, popsize=20, its = 1,000) and initial parameter sets defined on a log scale, and drawn from a uniform distribution on $[-6,6]^{3}$ for parameters $\beta$, $\alpha_{2}$ and $\gamma$, and on $[-1+\log(\hat{\mu}),\log(\hat{\mu})+1]$ for parameter $\mu$.

Confidence intervals

Pointwise likelihood-based confidence intervals are constructed thanks to the likelihood ratio test, following raue_structural_2009 ; for each estimated parameter $\hat{\theta}_{i}$, we select all the parameters $\theta_{i}=x$ such that:

[TABLE]

where $\Delta_{0.95}=\chi^{2}(0.95,1)=3.84$ is the $0.95$-quantile of the $\chi^{2}$ law with $1$ degree of freedom.

Model selection.

AIC and BIC analyses were performed to compare the submodels. The reader can refer to burnham_model_2003 (Chapter 6) for a detailed presentation of the rule of thumb, classically used to analyze the $\Delta^{AIC}_{i}:=AIC_{i}-AIC_{\min}$ and $\Delta^{BIC}_{i}=BIC_{i}-BIC_{\min}$ values, where $i$ is the index of the $i$th model:

• •

a $\Delta$ value lower than 2 indicates that the considered model is almost as probable as the “best” model;

• •

a $\Delta$ value between 2 and 7 suggests that the considered model is a suitable alternative to the “best” model;

• •

a $\Delta$ value between 7 and 10 suggests that the considered model is less relevant than the “best” model;

• •

a $\Delta$ value upper than 10 suggests that the considered model can be safely ruled out.

This $\Delta$ approach is completed by the AIC and BIC weight analyzes. For each dataset and criterion (AIC or BIC), we order the AIC/BIC weights from the highest to the lowest values. We then compute the cumulative sum of these weights, starting from the highest one. The selected models are the first ones such that the cumulative sum reaches the threshold p-value $0.95$.

6.7 Detailed calibration analysis

Two-event submodels

The fitting results obtained for submodels $(\mathcal{R}_{1},\mathcal{R}_{3})$ and $(\mathcal{R}_{1},\mathcal{R}_{4})$ from the experimental datasets are shown in Figure 5 and discussed in the main text, Section 4.3. One fitting result for the in silico datasets and for submodels $(\mathcal{R}_{1},\mathcal{R}_{3})$ and $(\mathcal{R}_{1},\mathcal{R}_{4})$ is shown in Figure 10. We verify that the inferred trajectories are coherent with the selected datasets.

In Figures 11, we show the PLE for each estimated parameter in each in-silico dataset. Both the initial condition parameter $\mu$ (orange solid lines) and asymmetric division rate $\alpha_{2}$ (green solid line) are practically identifiable (in the sense given in raue_structural_2009 ), while parameter $\gamma$ (blue solid line) is only partially practically identifiable in most cases. We observe that both parameters $\alpha_{2}$ ($\mathcal{R}_{3}$) and $\gamma$ ($\mathcal{R}_{4}$) are practically identifiable and close to their expected values (less than one $log10$ of difference) when the parameters are of the same order of magnitude than $\alpha_{1}$. In contrast, a small parameter value compared to $\alpha_{1}$ leads to a biased parameter estimate, with a huge shift between the estimated and true parameter values (up to two $log10$ difference).

The estimator for the initial condition parameter $\mu$ may also be slightly biased with submodel $(\mathcal{R}_{1},\mathcal{R}_{3})$ (less than one $log10$ of difference) compared to submodel $(\mathcal{R}_{1},\mathcal{R}_{4})$ .

Three-event submodels and complete model

We turn now to the analysis of three-event submodels $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3})$, $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{4})$ and $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$) and the complete model ($(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3},\mathcal{R}_{4})$. Qualitatively, the fitting results for submodel $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3})$ are similar to those for submodel $(\mathcal{R}_{1},\mathcal{R}_{3})$ (data not-shown); they are characterized by a high probability to produce ten or more proliferative cells before the precursor cell extinction. The fitting results for submodels $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{4})$ and $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$ are rather similar to submodel

$(\mathcal{R}_{1},\mathcal{R}_{4})$; they are characterized by direct cell transition with very little concomitant cell proliferation, followed by prolonged cell proliferation after precursor cell extinction. The fitting results for the complete model are shown in the bottom panels of Figure 5 for both the Wild-type and Mutant subsets and discussed in the main text, Section 4.3.

The PLEs for each dataset and each parameter are presented in Figure 12 for the three-event submodels. The corresponding parameter values and confidence intervals for the Wild-Type and Mutant subsets are given in Tables 4 and 5. As observed for the two-event submodels, in each case, the initial condition parameter $\mu$ (orange solid lines) is always practically identifiable, and its fitted value is close to the true one for the in silico datasets. In contrast, all other parameters have a lack of identifiability, both with the experimental and in silico datasets. Specifically, the asymmetric division rate $\alpha_{2}$ is practically not identifiable for submodel $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3})$ with the experimental subsets. Interestingly, when the asymmetric division event is combined with the symmetric division event (submodel $(\mathcal{R}_{1},\mathcal{R}_{3},\mathcal{R}_{4})$) rather than with the auto-amplified transition (submodel $(\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3})$), the asymmetric division rate $\alpha_{2}$ becomes identifiable in the experimental subsets, which reveals complex parameter dependencies between the asymmetric division rate $\alpha_{2}$ and auto-amplified transition rate $\beta$.

Acknowledgements.

The authors wish to thank Ken McNatty for providing the experimental dataset and Danielle Monniaux for helpful discussions.