Multiscale population dynamics in reproductive biology: singular perturbation reduction in deterministic and stochastic models

TL;DR

This paper explores multiscale modeling of ovarian follicle population dynamics using differential equations and singular perturbation techniques, revealing key population features and convergence behaviors.

Contribution

It introduces a unified approach applying singular perturbation reduction to both deterministic and stochastic models of follicle dynamics, highlighting population feedback effects.

Findings

Convergence of models to limit behaviors demonstrated.

Reproduction of bimodal population distribution.

Identification of a slope break in follicle decay with age.

Abstract

In this study, we describe different modeling approaches for ovarian follicle population dynamics, based on either ordinary (ODE), partial (PDE) or stochastic (SDE) differential equations, and accounting for interactions between follicles. We put a special focus on representing the population-level feedback exerted by growing ovarian follicles onto the activation of quiescent follicles. We take advantage of the timescale difference existing between the growth and activation processes to apply model reduction techniques in the framework of singular perturbations. We first study the linear versions of the models to derive theoretical results on the convergence to the limit models. In the nonlinear cases, we provide detailed numerical evidence of convergence to the limit behavior. We reproduce the main semi-quantitative features characterizing the ovarian follicle pool, namely a bimodal distribution of the whole population, and a slope break in the decay of the quiescent pool with aging.

Figures

Tables

Table 1. Table 1. Parameter values for the numerical simulations

$s$	$p_1$	$x_1$	$p_2$	$x_2$	$c$
0.1	0.7	0.2	0.3	0.7	2.61
$f_0$	$g_0$	$K_{1,0}$	$K_{2,0}$	${\rho }_{ini}^0$	$a$	$K_{1,d},d\ne 0$	$K_{2,d},d\ne 0$
1	0.1	2	0	100	$1_{[0.3,0.7]}$	0	0

Equations

$$\begin{array}{ccccccc}{F}_0& \stackrel{{\lambda }_0}{\to }& F_1& \stackrel{{\lambda }_1}{\to }& \cdots & \stackrel{{\lambda }_{N-1}}{\to }& F_N\\ \downarrow {\mu }_0& & \downarrow {\mu }_1& & \downarrow \cdots & & \downarrow {\mu }_N\\ \varnothing \varnothing \varnothing \varnothing & & & & & & \end{array}$$

$$\left\{\begin{array}[]{lll}\displaystyle\frac{dy_{0}(t)}{dt}&\displaystyle=&\displaystyle-\big{(}\bar{\lambda_{0}}(y(t))+\bar{\mu_{0}}(y(t))\big{)}y_{0}(t)\,,\\ \\ \displaystyle\frac{dy_{1}(t)}{dt}&\displaystyle=&\displaystyle\bar{\lambda_{0}}(y(t))y_{0}(t)-\big{(}\lambda_{1}(y(t))+\mu_{1}(y(t))\big{)}y_{1}(t),\\ \\ \displaystyle\frac{dy_{i}(t)}{dt}&\displaystyle=&\displaystyle\lambda_{i-1}(y(t))y_{i-1}(t)-\big{(}\lambda_{i}(y(t))+\mu_{i}(y(t))\big{)}y_{i}(t),\quad i\in\{2,\dots,d\}\end{array}\right.$$

$$\overline{{\lambda }_0}(y)=\frac{\overline{f_0}}{1+K_{1,0}\sum _{j=1}^d{a}_j{y}_j},\overline{{\mu }_0}(y)=\overline{g_0}\left(1+K_{2,0}\sum _{j=1}^d{b}_j{y}_j\right),$$

$${\lambda }_i(y)=\frac{f_i}{1+K_{1,i}\sum _{j=1}^d{\omega }_{1,j}{y}_j},{\mu }_i(y)=g_i\left(1+K_{2,i}\sum _{j=1}^d{\omega }_{2,j}{y}_j\right),$$

$$\frac{d}{dt}\sum _{i=0}^d{y}_i(t)=-\overline{{\mu }_0}(y(t))y_0(t)-\sum _{i=1}^d{\mu }_i(y(t))y_i(t)\le 0,$$

$$\overline{f_0}=\epsilon f_0,\overline{g_0}=\epsilon g_0,y_0^{in}=\frac{x_0^{in}}{\epsilon }$$

$$x_0(t)=\epsilon y_0(t/\epsilon ),\text{ and for all }i\ge 1,x_i(t)=y_i(t/\epsilon ).$$

$$\left\{\begin{array}[]{lll}\displaystyle\frac{dx_{0}(t)}{dt}&\displaystyle=&\displaystyle-\big{(}\lambda_{0}(x(t))+\mu_{0}(x(t))\big{)}x_{0}(t)\\ \\ \displaystyle\varepsilon\frac{dx_{i}(t)}{dt}&\displaystyle=&\displaystyle\lambda_{i-1}(x(t))x_{i-1}(t)-\big{(}\lambda_{i}(x(t))+\mu_{i}(x(t))\big{)}x_{i}(t),\quad i\in\{1,\dots,d\}\end{array}\right.$$

$${\lambda }_0(x)=\frac{f_0}{1+K_{1,0}\sum _{i=1}^d{a}_i{x}_i},{\mu }_0(x)=g_0\left(1+K_{2,0}\sum _{j=1}^d{b}_j{x}_j\right),$$

$$\frac{d}{dt}{x}_0(t)+\epsilon \frac{d}{dt}\sum _{i=1}^d{x}_i(t)=-\sum _{i=0}^d{\mu }_i(x(t))x_i(t)\le 0,$$

$$\left\{\begin{array}[]{lll}\displaystyle\frac{d\bar{x}_{0}(t)}{dt}&\displaystyle=&\displaystyle-\big{(}\lambda_{0}(\bar{x}(t))+\mu_{0}(\bar{x}(t))\big{)}\bar{x}_{0}(t)\\ \displaystyle 0&\displaystyle=&\displaystyle\lambda_{i-1}(\bar{x}(t))\bar{x}_{i-1}(t)-\big{(}\lambda_{i}(\bar{x}(t))+\mu_{i}(\bar{x}(t))\big{)}\bar{x}_{i}(t),\quad i\in\{1,\dots,d\}\end{array}\right.$$

$$\left\{\begin{array}[]{lll}\displaystyle\bar{x}_{0}(t)&\displaystyle=&\displaystyle x_{0}^{in}\exp\left(-(f_{0}+g_{0})t\right)\,,\\ \displaystyle\bar{x}_{i}(t)&\displaystyle=&\displaystyle\left(\prod_{j=0}^{i-1}\frac{f_{j}}{f_{j+1}+g_{j+1}}\right)\bar{x}_{0}(t),\quad i\in\{1,\dots,d\}\,.\end{array}\right.$$

$$\left\{\begin{array}[]{lll}\displaystyle\frac{d\bar{x}_{0}(t)}{dt}&\displaystyle=&\displaystyle-\big{(}\lambda_{0}(\bar{x}(t))+\mu_{0}(\bar{x}(t))\big{)}\bar{x}_{0}(t)\,,\\ \displaystyle\bar{x}_{1}(t)&\displaystyle=&\displaystyle\frac{-(f_{1}+g_{1})+\sqrt{(f_{1}+g_{1})^{2}+4f_{0}\bar{x}_{0}(t)(f_{1}+g_{1})K_{1,0}\sum_{i=1}^{d}\,a_{i}\prod_{j=1}^{i-1}\frac{f_{j}}{f_{j+1}+g_{j+1}}}}{2(f_{1}+g_{1})K_{1,0}\sum_{i=1}^{d}\,a_{i}\prod_{j=1}^{i-1}\frac{f_{j}}{f_{j+1}+g_{j+1}}}\,,\\ \displaystyle\bar{x}_{i}(t)&\displaystyle=&\displaystyle\left(\prod_{j=1}^{i-1}\frac{f_{j}}{f_{j+1}+g_{j+1}}\right)\bar{x}_{1}(t),\quad i\in\{2,\dots,d\}\,.\end{array}\right.$$

$$f_0{\bar{x}}_0=(f_1+g_1)\left(1+K_{1,0}\sum _{i=1}^d{a}_i\prod _{j=1}^{i-1}\frac{f_j}{f_{j+1}+g_{j+1}}{\bar{x}}_1\right){\bar{x}}_1$$

$$\left\{\begin{array}[]{lll}\displaystyle x_{0}^{\varepsilon}(t)&\displaystyle=&\displaystyle x_{0}^{in}\exp\left(-(f_{0}+g_{0})t\right),\\ \displaystyle(x_{1}^{\varepsilon},\cdots,x_{d}^{\varepsilon})^{T}(t)&\displaystyle=&\displaystyle e^{-\frac{B\,t}{\varepsilon}}(x_{1}^{in},\cdots,x_{d}^{in})^{T}+\int_{0}^{t}\frac{1}{\varepsilon}e^{-\frac{B(t-s)}{\varepsilon}}e_{1}f_{0}x_{0}^{\varepsilon}(s)ds,\end{array}\right.$$

$$B=\left(\begin{array}{cccc}(f_1+g_1)& 0& \dots & 0\\ -f_1& \ddots & & ⋮\\ 0& \ddots & \ddots & ⋮\\ ⋮& & & 0\\ 0& 0& -f_{d-1}& (f_d+g_d)\end{array}\right).$$

$$\underset{\epsilon \to 0}{\lim }\underset{t>\eta }{\sup }\underset{i\in \{1,\dots ,d\}}{\max }\mid x_i^{\epsilon }(t)-{\bar{x}}_i(t)\mid =0.$$

$${\int }_0^t\frac{1}{\epsilon }{e}^{-\frac{B(t-s)}{\epsilon }}{e}_1{f}_0{x}_0^{\epsilon }(s)ds={\int }_0^t{B}^{-1}{e}^{-\frac{B(t-s)}{\epsilon }}{e}_1{f}_0(f_0+g_0)x_0^{\epsilon }(s)ds+B^{-1}{e}_1{f}_0{x}_0^{\epsilon }(t)-B^{-1}{e}^{-\frac{Bt}{\epsilon }}{e}_1{f}_0{x}_0^{in}.$$

$$\begin{gathered} \Vert {\int }_0^t{B}^{-1}{e}^{-\frac{B(t-s)}{\epsilon }}{e}_1{f}_0(f_0+g_0)x_0^{\epsilon }(s)ds\Vert \le f_0(f_0+g_0)x_0^{in}\Vert {\int }_0^t{B}^{-1}{e}^{-\frac{B(t-s)}{\epsilon }}{e}_{1}ds\Vert \\ \le f_0(f_0+g_0)x_0^{in}\epsilon \Vert B^{-2}(Id-e^{-\frac{Bt}{\epsilon }})e_1\Vert \le f_0(f_0+g_0)x_0^{in}\epsilon |||B^{-2}|||(1+e^{-\min (f_i+g_i)\frac{t}{\epsilon }}) \end{gathered}$$

$$\forall t\ge 0,\Vert e^{-\frac{Bt}{\epsilon }}{e}_1\Vert \le \Vert e_1\Vert e^{-\min (f_i+g_i)\frac{t}{\epsilon }}$$

$$\begin{gathered} \underset{t>\eta }{\sup }\Vert (x_1^{\epsilon },\cdots ,x_d^{\epsilon })(t)-({\bar{x}}_1,\cdots ,{\bar{x}}_d)(t)\Vert \le \left(\Vert (x_1^{in},\cdots ,x_d^{in})\Vert +|||B^{-1}|||f_0{x}_0^{in}\right)e^{-\min (f_i+g_i)\frac{\eta }{\epsilon }} \\ +2f_0(f_0+g_0)x_0^{in}|||B^{-2}|||\epsilon , \end{gathered}$$

$$E_1(t)=\sum _{i=0}^d\mid X_i^{\epsilon }(t)-{\bar{X}}_i(t)\mid ,$$

$$E_2(T)={\int }_0^T\sum _{i=0}^d\mid X_i^{\epsilon }(t)-{\bar{X}}_i(t)\mid dt,$$

$${\tilde{E}}_2(T)=\sum _{k=0}^{N_T}{\delta }_t\sum _{i=0}^d\mid X_i^{\epsilon }(t_k)-{\bar{X}}_i(t_k)\mid ,$$

$$\left\{\begin{array}[]{lll}\displaystyle\dfrac{d\rho_{0}(t)}{dt}&\displaystyle=&\displaystyle-(\lambda_{0}(\rho(t,.))+\mu_{0}(\rho(t,.)))\rho_{0}(t)\,,\\ \displaystyle\varepsilon\partial_{t}\rho(t,x)&\displaystyle=&\displaystyle-\partial_{x}(\lambda(\rho(t,.),x)\rho(t,x))-\mu(\rho(t,.),x)\rho(t,x)\,,\quad\mbox{for }x\in(0,1)\,,\\ \displaystyle\lim_{x\rightarrow 0}\lambda(\rho(t,.),x)\rho(t,x)&\displaystyle=&\displaystyle\lambda_{0}(\rho(t,.))\rho_{0}(t)\,,\\ \end{array}\right.$$

$${\lambda }_0(\rho (t,.))=\frac{f_0}{1+K_{1,0}{\int }_0^{1}a(y)\rho (t,y)dy},{\mu }_0(\rho (t,.))=g_0\left(1+K_{2,0}{\int }_0^{1}b(y)\rho (t,y)dy\right),$$

$$\lambda (\rho (t,.),x)=\frac{f(x)}{1+K_1(x){\int }_0^1{\omega }_1(y)\rho (t,y)dy},\mu (\rho (t,.),x)=g(x)\left(1+K_2(x){\int }_0^1{\omega }_2(y)\rho (t,y)dy\right)$$

$${\rho }_0(t=0)={\rho }_0^{in},\rho (t=0,x)={\rho }^{in}(x),x\in (0,1).$$

$$\frac{d}{dt}{\rho }_0(t)+\epsilon \frac{d}{dt}{\int }_0^1\rho (t,x)dx=-{\mu }_0(\rho (t,.)){\rho }_0(t)-{\int }_0^1\mu (\rho (t,.),x)\rho (t,x)dx-\underset{x\to 1}{\lim }\lambda (\rho (t,.),x)\rho (t,x)\le 0.$$

$$\left\{\begin{array}[]{lll}\displaystyle\partial_{t}\bar{\rho}_{0}(t)&\displaystyle=&\displaystyle-(\lambda_{0}(\bar{\rho}(t,.))+\mu_{0}(\bar{\rho}(t,.)))\bar{\rho}_{0}(t)\,,\\ \displaystyle\partial_{x}\big{(}\lambda(\bar{\rho}(t,.),x)\bar{\rho}(t,x)\big{)}&\displaystyle=&\displaystyle-\mu(\bar{\rho}(t,.))\bar{\rho}(t,x)\,,\quad\mbox{for }x\in(0,1)\,,\\ \displaystyle\lim_{x\rightarrow 0}\lambda(\bar{\rho}(t,.),x)\bar{\rho}(t,x)&\displaystyle=&\displaystyle\lambda_{0}(\bar{\rho}(t,.))\bar{\rho}_{0}(t)\,,\end{array}\right.$$

Full text

{resume}

Dans cette étude, nous décrivons différentes approches de modélisation de la dynamique des populations de follicules ovariens, basées sur des équations différentielles ordinaires (EDO), aux dérivées partielles (EDP) ou stochastiques (SDE), et tenant compte des interactions entre follicules. Nous avons mis un accent particulier sur la représentation des rétro-actions exercées par les follicules en croissance sur l’activation des follicules quiescents. Nous tirons parti de la différence d’échelle de temps entre les processus de croissance et d’activation pour appliquer des techniques de réduction de modèle dans le cadre des perturbations singulières. Nous étudions d’abord les versions linéaires du modèle afin d’en déduire des résultats théoriques sur la convergence vers le modèle limite. Dans le cas non linéaire, nous fournissons des arguments numériques détaillés sur la convergence vers le comportement limite. Nous reproduisons les principales caractéristiques semi-quantitatives caractérisant le pool de follicules ovariens, à savoir une distribution bimodale de la population totale et une rupture de pente dans la décroissance du pool de follicules quiescents avec le vieillissement.

Introduction

In mammals, the pool of oocytes (egg cells) available for a female throughout her reproductive life is fixed very early, either during the fetal life or in the perinatal period. All along their maturation, oocytes are sheltered within spheroidal somatic structures called ovarian follicles. Folliculogenesis is the process of growth and maturation undergone by ovarian follicles from the time they leave the pool of quiescent, primordial follicles until ovulation, when they release a fertilizable oocyte.

Follicle growth is first due to the enlargement of the oocyte, then to the proliferation of somatic cells organized into successive concentric cell layers, and finally to the inflation of a fluid-filled cavity (antrum) that forms above a critical size. The activation of primordial follicles can occur at any time once they are formed [15], even if they can remain quiescent for up to tens of years [19]. Growing follicles can progress along the first developmental stages (known as “basal development”) before puberty. The final developmental stages (known as “terminal development”) can only occur after puberty ; they are related to the dynamics of ovarian cycles, involving endocrine feedback loops between the ovaries on one side, and the hypothalamus and pituitary gland on the other side. The whole sequence of development spans several months, as assessed by cell kinetics studies [10] or grafts of ovarian cortex [3]. The terminal stages are the shortest ; they cover a few weeks at most.

Since follicle activation is asynchronous, all developmental stages can be observed in the ovaries at any time during reproductive life. The follicle distribution (mostly studied using the size as a maturity marker) has a characteristic bimodal pattern, which is remarkably preserved between species. This pattern remains similar with ovarian aging, yet with a decreased amplitude [4], as a result of the progressive exhaustion of the quiescent pool. Such a distribution is shaped not only by the differences in the follicle activation times, but also by the hormonal interactions between follicles [17]. In particular, the activation and growth rates in the earliest stages are moderated by the Anti-Müllerian Hormone (AMH) secreted locally by the subpopulation of “intermediary” follicles (rigorously speaking: from the fully activated one-layer stage to the pre-antral and small antral stages) [20]. At the other end, the selection of ovulatory follicles results from a competition-like process operating amongst terminally developing follicles [1], which is mediated by endocrine controls and associated with a species-specific number of ovulations. Namely, inhibin (a peptid hormone) and estradiol, produced by the mature follicles, feedback onto the pituitary gland, leading to a drop in a pituitary hormone (the Follicle-Stimulating Hormone) supporting follicle survival.

Less than one in a thousand of the follicles manage to reach the ovulatory stage. All others disappear through a degeneration process (atresia) associated with the death of the somatic cells (during mid and terminal folliculogenesis) or oocyte (in the quiescent pool and during early folliculogenesis). For instance, in humans, the quiescent pool size is of the order of 1 million follicles, amongst which only some hundreds will reach ovulation [16].

Experimentalist investigators have proposed classifying follicle development into different stages, according to morphological and functional criteria such as follicle and oocyte diameters, number of cell layers, number of somatic cells, antrum formation [10, 18]. Hence, a natural formalism to consider when modeling follicle population dynamics is that of compartmental modeling, using either deterministic or stochastic rates for transfer ($\lambda_{i}$) and exit ($\mu_{i}$) rates (see Equation (1)). Pioneering studies (see e.g. [7]) have focused on fitting the parameters entering these rates according to follicle numbers available in each developmental stage. However, these studies remained rather descriptive and considered at best possibly time-varying (piecewise constant) rates [8], yet with no follicle interaction.

[TABLE]

Most of the classification criteria change in a continuous manner. In addition, the most common variable available to monitor follicle development on the ovarian scale, follicle size, is an intrinsically continuous variable. Hence, another suitable modeling formalism is that of PDE models for structured population dynamics. Although the interest of such a formalism has been pointed out quite early [14], it has yet not been implemented.

Finally, in some situations, a discrete stochastic formalism can be useful both to handle finite-size effects and follow individual follicle trajectories. This is especially true for the relatively small cohort of terminally developing follicles, and for transient physiological regimes when follicle pools are either still replenishing, or, on the contrary getting progressively exhausted. In any case, such a formalism gives insight into the fluctuations around the average deterministic behavior.

In this study, we describe different modeling approaches for follicle population dynamics, based on either ODE, PDE or SDE, and accounting for interactions between follicles. We put a special focus on representing the population-level modulation exerted by growing ovarian follicles on the activation of quiescent, primordial follicles. We take advantage of the timescale difference between the growth and activation processes to apply model reduction techniques in the framework of singular perturbations (slow/fast systems).

The paper is organized as follows. We successively introduce the ODE, PDE and SDE formulation of the model for follicle population dynamics. We describe the initial (non-rescaled) model in the ODE case. In each case, we introduce (i) the model in rescaled timescale exhibiting a slow/fast structure with a small perturbation parameter ($\varepsilon>0$) and (ii) the model in the limit $\varepsilon\to 0$. We discuss the well-posedness of the limit models in two situations: the linear formulation and a weakly nonlinear formulation in which only the quiescent follicle population is subject to a feedback from the remaining of the population. In the linear case, we prove the convergence of the rescaled to the limit models. In the nonlinear case, we provide detailed numerical evidence of convergence. The numerical illustrations are settled within a biologically-realistic framework, allowing us to reproduce the main semi-quantitative features characterizing the dynamics of the ovarian follicle pool, namely a bimodal distribution of the whole population and a slope break in the decay of the quiescent pool with aging.

1. Compartmental, ODE-based model

1.1. Initial model

Starting from the schematic model (see Eq. (1)), we formalize a system of nonlinear ordinary differential equations (ODE) as follows. Let $d\in{\mathbb{N}}^{*}$ and $y=(y_{0},\dots,y_{d})$ be a function such that, for all $i\in\{0,\dots,d\}$, $y_{i}:t\in\mathbb{R}_{+}\mapsto y_{i}(t)$ represents the time evolution of the number of follicles of maturity $i$. Follicles in the first compartment ($i=0$) are named quiescent follicles, and their maturation and death rates are denoted by $\bar{\lambda_{0}}$ and $\bar{\mu_{0}}$, respectively. Follicles in the intermediate compartments ($1\leq i\leq d-1$) are the growing follicles, and may either mature and go to the next maturation stage $i+1$, at rate $\lambda_{i}$, or die at rate $\mu_{i}$. Follicles in the last compartment ($i=d$) are named the mature follicles and can only die at rate $\mu_{d}$, i.e. $\lambda_{d}=0$ (death in this compartment corresponds to either degeneration or ovulation). All the rates $(\mu_{i},\lambda_{i})$ may depend on the growing and mature follicles population (non-local interactions), which leads to the following nonlinear ODE system

[TABLE]

where, for $i=0$,

[TABLE]

with non-negative parameter constants $\bar{f_{0}}$, $\bar{g_{0}}$, $K_{1,0}$, $K_{2,0}$, and $a_{j}\in[0,1],b_{j}\in[0,1]$.

For $i\in\{1,\dots,d\}$,

[TABLE]

with non-negative parameter constants, $f_{i}$ ($f_{d}=0$), $K_{1,i}$, $g_{i}$, $K_{2,i}$ and $\omega_{1,j}\in[0,1]$, $\omega_{2,j}\in[0,1]$. The specific functional forms of the rate coefficients are motivated by biological knowledge (see Introduction). The population feedback tends to lower the maturation rate and to raise the death rate. With a non-negative vector $y^{in}\in{\mathbb{R}}_{+}^{d+1}$ as initial data, one can see that Eq. (2) generates a unique non-negative solution for all times (the right-hand side is globally Lipschitz, with positive off-diagonal entries). Moreover, one can obtain immediately the following conservation law,

[TABLE]

which shows that any follicle sub-population $y_{i}$ is globally bounded.

1.2. Rescaled model

As outlined in the Introduction, before reproductive senescence, quiescent follicles are very numerous compared to the growing and mature follicles, follicle activation dynamics are much slower than growth dynamics, yet the flow of follicles between each compartment is of the same order. In consistency with this timescale contrast, we introduce a small positive parameter $\varepsilon\ll 1$, such that

[TABLE]

with non-negative constants $f_{0}$, $g_{0}$ and positive initial data $x_{0}^{in}$, independent of $\varepsilon$. Note that the initial flow $\bar{f_{0}}y_{0}^{in}=f_{0}x_{0}^{in}$ is preserved.

We then define the rescaled solution $x=(x_{0},\dots,x_{d})$ by, for all $t\geq 0$,

[TABLE]

Then $x$ is solution of the following system

[TABLE]

with initial condition given by $x_{0}(t=0)=x_{0}^{in}$, and $x_{i}(t=0)=x_{i}^{in}:=y_{i}^{in}$, for $i\in\{1,\dots,d\}$, and where, for $i=0$,

[TABLE]

and $\lambda_{i}$ and $\mu_{i}$, for $i\in\{1,\dots,d\}$, are defined in Eq. (4). We note that the conservation law (5) becomes

[TABLE]

We now consider the limit for which the small parameter $\varepsilon$ tends to [math] and the associated sequence $(x^{\varepsilon})$ solution of system (8). In such a case, system (8) is called a “slow-fast” system ($x_{0}^{\varepsilon}$ is the slow variable, $(x_{1}^{\varepsilon},\dots,x_{d}^{\varepsilon})$ are the fast variables) and the study of the limit behavior when $\varepsilon\to 0$ is a singular perturbation problem (see for instance [21]).

1.3. Limit model

Formally, setting $\varepsilon=0$ in system (8) leads to the following system:

[TABLE]

with initial condition given by $\bar{x}_{0}(t=0)=x_{0}^{in}$, and $\bar{x}_{i}(t=0)=\bar{x}_{i}^{in}\geq 0$, that satisfies the second line of Eq. (11) at $t=0$. Note that system (11) is a differential-algebraic system, in which the variable $(\bar{x}_{1},\dots,\bar{x}_{d})$ can be seen as reaching instantaneously (at any time $t$) a quasi-steady state, “driven” by the time-dependent variable $\bar{x}_{0}(t)$.

System (11) is not necessarily well-posed, as there may be several solutions to the second line of Eq. (11). In the next two specific examples, we can prove that system (11) does admit a single positive solution, which is a natural limit candidate for the sequence $x^{\varepsilon}$.

Example 1** (Linear case).**

Let us suppose that $K_{1,i}=K_{2,i}=0$ for all $i\in\{0,\dots,d\}$, and $f_{i}+g_{i}>0$ for all $i\in\{1,\dots,d\}$. Then, system (11) becomes linear, and has a unique solution given by $\bar{x}=(\bar{x}_{0},\dots,\bar{x}_{d})$ such that for all $t\geq 0$,

[TABLE]

Example 2** (Feedback onto quiescent follicle activation and death rates).**

Let us suppose that $K_{1,i}=0$, $K_{2,i}=0$ and $f_{i}+g_{i}>0$ for $i\in\{1,\dots,d\}$, yet $K_{1,0}>0$ and $K_{2,0}\geq 0$. Then, system (11) with positivity requirement ($x_{i}\geq 0$ for $i\in\{0,\dots,d\}$) can be rewritten as:

[TABLE]

which admits a unique solution. Indeed, one can verify that $\bar{x}_{1}$ is the only positive root of a polynomial of degree $2$, namely

[TABLE]

1.4. Convergence in the linear case

In this paragraph, we assume that $K_{1,i}=K_{2,i}=0$ for all $i\in\{0,\dots,d\}$, and $f_{i}+g_{i}>0$ for all $i\in\{1,\dots,d\}$ as in Example 1. In such a case, one can solve system (8) explicitly for each $\varepsilon>0$. The solution is given by, using vectorial notations, for all $t\geq 0$,

[TABLE]

where $e_{1}=(1,0,\cdots,0)^{T}$ and

[TABLE]

It is thus clear that $x_{0}^{\varepsilon}=\bar{x}_{0}$ is a constant sequence in $\varepsilon$ (as both $x_{0}^{\varepsilon}$ and $\bar{x}_{0}$ have same initial conditions, and same evolution equation). For the fast variables, we prove the following

Proposition 1**.**

Assume that $K_{1,i}=K_{2,i}=0$ for all $i\in\{0,\dots,d\}$, and $f_{i}+g_{i}>0$ for $i\in\{1,\dots,d\}$. Then, for all $\eta>0$, we have

[TABLE]

where $(\bar{x}_{1},\cdots,\bar{x}_{d})$ is given in Eq. (12).

Proof.

From Eq. (14) and initial condition, and using integration by parts, we obtain

[TABLE]

As $x_{0}^{\varepsilon}$ is uniformly bounded (both in $\varepsilon$ and time) by $x_{0}^{in}$, we obtain, taking the $1$-norm,

[TABLE]

with $|||B^{-2}|||=\displaystyle\sup_{x\in\mathbb{R}^{d},x\neq 0}\frac{\left\lVert B^{-2}x\right\rVert}{\left\lVert x\right\rVert}$. The third inequality above was deduced from

[TABLE]

We verify that for all $t\geq 0$, $B^{-1}e_{1}f_{0}x_{0}^{\varepsilon}(t)=(\bar{x}_{1},\cdots,\bar{x}_{d})(t)$. Then, we obtain,

[TABLE]

which converges to [math] as $\varepsilon$ converges to [math]. ∎

Remark 1**.**

The proof of Proposition 1 can also be obtained as a direct application of Tikhonov theorem [21].

Remark 2**.**

It is apparent in formula (14) that one cannot hope to obtain a convergence on a time interval starting from [math] (unless the initial data is “well-prepared”), and that standard Ascoli-Arzela theorem would not apply in such a case, as the time derivative of $(x_{1}^{\varepsilon},\cdots,x_{d}^{\varepsilon})(t)$ is not uniformly bounded as $\varepsilon\to 0$.

1.5. Numerical convergence

In this paragraph, we illustrate the convergence of $(x_{0}^{\varepsilon},x_{1}^{\varepsilon},\cdots,x_{d}^{\varepsilon})$ towards $(\bar{x}_{0},\bar{x}_{1},\cdots,\bar{x}_{d})(t)$ in a nonlinear scenario. The chosen scenario and the parameter values are detailed in the Appendix (section 4).

In Figure 1, we plot the trajectories $x_{i}(t)$ in each maturity compartment ($d=10$) for the rescaled and limit models on a time horizon $t\in(0,1)$. In each compartment, the trajectories of the rescaled model get closer and closer to the limit model as $\varepsilon\to 0$. For $\varepsilon=0.001$, they are almost indistinguishable. Note however that, for the growing follicles, the initial conditions of the rescaled and limit models are different, and the convergence holds only for positive times.

In Figure 2, using the same parameters as in Figure 1, we display the maturity distribution in the growing follicle population, for various $\varepsilon$. We can expect from Figure 1 that the convergence gets better for larger times. We confirm this fact in Figure 2 where we compare the maturity distribution in the growing follicle population at two times, $t=0.1$ and $t=1$. We can quantify the error between the rescaled model and the limit model by computing the $l^{1}$ error at time t,

[TABLE]

and the cumulative error on time interval $[0,T]$,

[TABLE]

which can be assessed numerically as

[TABLE]

where $t_{k}=k\delta_{t}$, for k=0{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\cdots}N_{T}. From the right panel of Figure 2, we can see that, at least for a small enough $\varepsilon$, the error is inversely proportional to $\varepsilon$, $E_{i}\approx C^{te}\varepsilon^{-1}$, where the constant pre-factor may depend on the chosen norm or the particular time $t$.

2. PDE model

When considering a continuous maturity variable, the PDE formalism is more suited for representing the follicle population dynamics. In this section, we skip the rescaling procedure, which follows an analogous reasoning as that detailed in section 1, and present directly the rescaled model.

2.1. Rescaled model

Denoting by $\rho_{0}(t)$ the number of quiescent follicles and by $\rho(t,x)$ the population density of follicles of maturity $x$, we consider the following coupled ODE-PDE system, for all $t\geq 0$,

[TABLE]

where

[TABLE]

and, for all $x\in(0,1)$,

[TABLE]

with initial condition

[TABLE]

We assume that $f,g,K_{1},K_{2},a,b,w_{1},w_{2},\rho^{in}$ are regular enough functions, and will admit the existence and uniqueness of solutions of system (20). A standard fixed point argument, based on the mild formulation, could be used (see for instance [2, 6, 5] and references therein), yet this is beyond the scope of this work. We can write the following conservation law, that gives (at least formally)

[TABLE]

In the following, we consider a sequence $(\rho_{0}^{\varepsilon},\rho^{\varepsilon})$ of solutions of system (20) in the limit $\varepsilon\to 0$.

2.2. Limit model

Formally, setting $\varepsilon=0$ in system (20) leads to the following system: for all $t\geq 0$,

[TABLE]

with an initial condition given by $\bar{\rho}_{0}(t=0)=\rho_{0}^{in}$, and $\bar{\rho}(t=0,.)=\bar{\rho}^{in}$, that satisfies the second and third lines of Eq. (25) at $t=0$. System (25) is not necessarily well-posed, as there may be several solutions $\bar{\rho}$ for a given $\bar{\rho}_{0}$. In the next two specific examples, we can prove that system (25) does admit a single positive solution, which is a natural limit candidate for the sequence $(\rho_{0}^{\varepsilon},\rho^{\varepsilon})$.

Example 3** (Linear case).**

Let us suppose that $K_{1,0}=K_{2,0}=0$ and $K_{1}\equiv 0$, $K_{2}\equiv 0$. Assume furthermore that $f(0)>0$. Then, system (25) becomes linear, and has a unique solution given by

[TABLE]

Example 4** (Feedback onto quiescent follicle activation and death rates).**

Let us suppose that $K_{1}\equiv 0$, $K_{2}\equiv 0$, and $K_{1,0}>0$ and $K_{2,0}\geq 0$. Assume furthermore that $f(0)>0$. Then, system (25) with positivity requirement $\bar{\rho}(t,x)\geq 0$ can be rewritten as

[TABLE]

which admits a unique solution. Indeed, the functional expression of $\bar{\rho}$ comes directly from solving the second line of Eq. (25). Using the boundary condition in the third line of Eq. (25), one can verify that $\bar{\rho}(t,0)$ is the positive root of a polynomial of degree $2$, namely

[TABLE]

The limit system (27) is a nonlinear ODE. To simplify notations, we introduce

[TABLE]

from which we rewrite system (27) as

[TABLE]

We can thus solve $\bar{\rho}_{0}(t)$ as the solution of an autonomous nonlinear ODE

[TABLE]

with

[TABLE]

from which we then compute $\bar{\rho}(t,0)$ and eventually $\bar{\rho}(t,x)$.

2.3. Convergence in the linear case

In this paragraph, we assume that $K_{1,0}=K_{2,0}=0$, $K_{1}\equiv 0$, $K_{2}\equiv 0$, and $f(0)>0$ as in Example 3. In such a case, one can solve explicitly system (20) for each $\varepsilon>0$ using the characteristics method. We obtain

[TABLE]

where $X(0;t,x)$ is the location of the characteristic at time [math], given that it goes through the point $x$ at time $t$, namely:

[TABLE]

It is thus clear that $\rho_{0}^{\varepsilon}=\bar{\rho}_{0}$ is a constant sequence in $\varepsilon$. For the population density $\rho^{\varepsilon}$, which is here the fast unknown, we prove the following

Proposition 2**.**

Assume that $K_{1,0}=K_{2,0}=0$, $K_{1}\equiv 0$, $K_{2}\equiv 0$, and $f(0)>0$ with $\int_{0}^{1}\frac{1}{f(x)}dx<\infty$ and $\int^{1}_{0}\frac{g(x)+f^{\prime}(x)}{f(x)}dx<\infty$. Then, for all $\eta>0$, we have

[TABLE]

where $\bar{\rho}$ is given in Eq. (26).

Proof.

It is clear that for any $\eta>0$, there exists $\varepsilon^{\prime}$ such that for all $\varepsilon<\varepsilon^{\prime}$, and all $t>\eta$ we have

[TABLE]

Then, comparing the solutions of Eq. (29) and Eq. (26) allows us to conclude that, for all $t>\eta$,

[TABLE]

∎

Remark 3**.**

As in Remark 1, we can see that during a time of order $\varepsilon t$, we cannot get the convergence of the rescaled model towards the reduced one, which precludes uniform convergence in time starting from {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}t=}0.

2.4. Numerical study

In this paragraph, we detail the numerical schemes that we have designed to solve both systems (20) and (25), and illustrate the consistency and convergence of these algorithms using the exact solutions.

2.4.1. Numerical scheme for the limit model

We design a finite difference scheme to compute a numerical solution to the PDE limit system (25). This system is nonlinear due to the dependence of $\lambda_{0}$, $\lambda$ and $\mu$ upon the solution $\bar{\rho}(t,x)$, which itself depends on $\bar{\rho}_{0}$. We propose to treat this nonlinearity with a fixed point scheme. At each time step $t_{n}=n\Delta_{t}$, for $n=0,\ldots,N$ with $T=N\Delta_{t}$ we build a sequence $\rho^{\ell}(t_{n},x)$ such that

[TABLE]

Let $x_{k}=k\Delta_{x}$, for $k=0,\ldots,M$, with $M\Delta_{x}=1$. We introduce the discretized approximations

[TABLE]

which we compute as follows. Let $\eta\ll 1$.

1. (1)

Initialization.

[TABLE] 2. (2)

For $n=0\nearrow N$ compute $\bar{\rho}_{k}^{n}=\lim_{\ell\rightarrow\infty}\rho_{k}^{\ell,n}$ iteratively then update $\bar{\bar{\rho}}_{0}^{n+1}$ as follows.

1. (a)

Initialize residual $R_{\ell}^{n}=1$ and set $\ell=0$ 2. (b)

While $R_{\ell}^{n}>\eta$ do

• •

Compute the PDE parameters $\bar{\lambda}_{0}^{n,\ell}$, $\lambda_{k}^{n,\ell}$ and $\mu_{k}^{n,\ell}$ (for $k=0,\ldots,M$) by standard trapezoidal rules

[TABLE]

• •

Enforce boundary condition at $x=0$

[TABLE]

• •

Integrate numerically the PDE in $x$

[TABLE]

• •

Compute residual between $\ell$ and $\ell+1$ iterations

[TABLE] 3. (c)

Set the new population density equal to the value obtained at the end of the fixed point procedure $\bar{\rho}_{k}^{n}=\rho_{k}^{n,\ell+1}$ 4. (d)

Compute the ODE parameters from the fixed point value, using standard trapezoidal rules

[TABLE] 5. (e)

Integrate numerically the ODE $\bar{\bar{\rho}}_{0}$ between $t_{n}$ and $t_{n+1}$ with a classic explicit Euler scheme

[TABLE]

2.4.2. Convergence of the numerical scheme for the limit model.

For the nonlinear scenario described in Example 4 we obtain a “pseudo exact” solution using scipy library ODE solver odeint to solve Eq. (28) and compute $\bar{\rho}_{0}(t)$, and we use Eq. (27) to compute $\bar{\rho}(t,\cdot)$. We use this reference solution to assess the performances of the numerical scheme described in paragraph 2.4.1. We display on Figure 3 the simulations performed with the parameter values detailed in the Appendix (section 4), on a time horizon $t\in(0,1)$, and an increasingly fine discretization

[TABLE]

The top left panel shows the difference between the pseudo exact solution $\bar{\rho}_{0}(t_{n})$ of Eq. (28) and the numerical solutions $\bar{\bar{\rho}}^{n}$ computed with the fixed point algorithm, as a function of $t\in[0,T]$. The bottom left panel shows the pseudo exact solution $\bar{\rho}(T,x_{k})$ and the numerical solutions $\bar{\rho}^{N}_{k}$ at final time $T=1$, as a function of $x\in[0,1]$. In the top right panel, we display the relative errors in $L^{1}$ and $L^{\infty}$ norms between the pseudo exact and the numerical solution $\bar{\rho}_{0}(t)$, as a function of $\Delta_{t}=\Delta_{x}=1/N$:

[TABLE]

[TABLE]

and the relative errors between the pseudo exact and the numerical solutions $\bar{\rho}(T,x)$

[TABLE]

[TABLE]

As expected, the errors are linear in $\Delta_{x}$, which means that the order of numerical convergence is one. In the bottom right panel, we display the number of iterations performed in the inner loop of the fixed point algorithm, as a function of time. The number of iterations needed to converge decreases with time and the number of time points. This tendency is not really surprising, since, at each time step, we start with the solution obtained at the previous time step as initial condition for the fixed point loop. Since the solution decreases with time, the distance between the fixed point initial condition and the solution decreases with both the time and time step, hence convergence requires less iterations.

2.4.3. Numerical scheme for the rescaled model

We design an explicit finite volume scheme to compute a numerical solution to the rescaled model (Eq. 20). The discretized unknowns are at each time step $t_{n}=n\Delta_{t}$, for $n=0,\ldots,N$ with $T=N\Delta_{t}$

[TABLE]

and, for $x_{k}=k\Delta_{x}$, $k=0,\ldots,M$, with $M\Delta_{x}=1$,

[TABLE]

We integrate numerically the PDE between $t_{n}$ and $t_{n+1}$ and over $[{x_{k}},{x_{k+1}}]$ by freezing the nonlinear coefficients $\lambda(\bar{\rho}(t,.),x)$ and $\mu(\bar{\rho}(t,.))$ at time $t_{n}$

[TABLE]

At each time step, we compute both the PDE and ODE coefficients using the midpoint rule and the numerical solution $(\rho^{\varepsilon,n}_{k})_{k=0,\ldots,M}$ as a piecewise constant solution

[TABLE]

For the explicit scheme (32), two stability conditions must be satisfied

• •

CFL-like stability condition:

[TABLE]

which can be rewritten as

[TABLE]

• •

Positivity conservation condition: $1+\frac{\Delta_{t}^{n}}{\varepsilon}\left(\frac{\lambda_{k}^{\varepsilon,n}}{\Delta_{x}}-\mu_{k}^{\varepsilon,n}\right)\geq 0,\quad k=1,\ldots,M$ we impose that if $\rho_{i}^{\varepsilon,n}=\delta_{ik}$ then $\rho_{i}^{\varepsilon,n+1}\geq 0$ for all $k$ and all $i$ which leads to

[TABLE]

which can be rewritten as

[TABLE]

The overall numerical scheme proceeds as follows:

1. (1)

Initialization :

[TABLE] 2. (2)

For $n=0\nearrow N$ compute $\rho_{k}^{\varepsilon,n+1}$ for $k=1,\ldots,M$ then update $\bar{\rho}_{0}^{\varepsilon,n+1}$ and finally compute $\rho_{0}^{\varepsilon,n+1}$ :

1. (a)

Compute the PDE and ODE coefficients $\bar{\lambda}_{0}^{\varepsilon,n}$, $\bar{\mu}_{0}^{\varepsilon,n}$, $\lambda_{k}^{\varepsilon,n}$ and $\mu_{k}^{\varepsilon,n}$ (Eq. 37) 2. (b)

Compute $\Delta_{t}^{n}$ satisfying stability conditions (38) and (39) 3. (c)

Integrate numerically the PDE in $x$ at time $t_{n}$ using (32) 4. (d)

Integrate numerically the ODE ${\bar{\rho}}_{0}$ between $t_{n}$ and $t_{n+1}$

[TABLE] 5. (e)

Enforce the boundary condition at $x=0$

[TABLE]

2.4.4. Convergence of the numerical scheme for the rescaled model.

We start by checking the convergence of the numerical scheme in a case where we know an exact solution, that is the linear case $K_{10}=K_{20}=K_{1}=K_{2}=0$. We test several discretizations

[TABLE]

for the parameter values of the linear scenario detailed in the Appendix (section 4). The results are displayed in Figure 4. In the left panel, we see that $\rho_{0}(t)$ is computed exactly since all curves corresponding to different discretizations are superimposed (as expected in the linear case). In the center panel, we display the solution at a fixed time $T$ as a function of $x$, which does depend on the discretization. The right panel shows the relative error curves, which exhibit a convergence rate better than linear.

2.4.5. $\varepsilon$-convergence towards the limit model

So far we have proved the convergence of the rescaled model towards the limit model when $\varepsilon\to 0$ in the linear case. We can only test it numerically in the general case. To minimize the numerical error arising from solving the limit model numerically, we illustrate the $\varepsilon$-convergence of $(\rho_{0}^{\varepsilon},\rho^{\varepsilon})$ towards $(\bar{\rho}_{0},\bar{\rho})$ in the nonlinear scenario of example 4 (see details on the parameter values in the Appendix (section 4)), for which the pseudo-exact solution of the limit model is available. As in Figure 3, pseudo exact solutions $(\bar{\rho}_{0},\bar{\rho})$ are simulated using Eq. (28) and Eq. (27). The results are displayed in Figure 5 for a set of $\varepsilon$ values and two discretizations $M=N=100$ and $M=N=200$. The agreement of $\rho_{0}(t)$ with the limit model solution (top left panel), and that of $\rho^{\varepsilon}(t,x)$ (bottom left panel), are qualitatively good as soon as $\varepsilon\leq 10^{-2}$. The relative error for $\rho_{0}(t)$ (top right panel) exhibits a linear behavior in $\varepsilon$. The error curves for the convergence of the solution $\rho^{\varepsilon}(t,x)$ in the domain (bottom right panel) are not linear, and remain beyond a threshold when $\varepsilon$ goes to 0. However, both the value of $\varepsilon$ for which the error approaches this threshold, and the value of the error itself decrease when we refine the discretization. This indicates that we should refine the discretization when we decrease $\varepsilon$. Since our current numerical scheme is explicit in time, any refinement must be simultaneous in $\Delta_{x}$ and $\Delta_{t}$, and as a consequence the cost in CPU time depends quadratically in $\varepsilon^{-1}$.

In practice, realistic values for $\varepsilon$ should remain tractable. However this behavior is a good incentive to study a more economical numerical scheme, namely an implicit one, which would provide accurate results with coarser discretizations.

3. SDE model

We now turn to a stochastic model for the follicle population dynamics. We skip the rescaling procedure, which follows an analogous reasoning as that detailed in section 1, and present directly the rescaled model.

3.1. Rescaled model

We consider the following coupled Poisson-driven SDE system, given by, for all $t\geq 0$,

[TABLE]

where $X^{\varepsilon}$ represents the vectorial process $(X_{0}^{\varepsilon},\cdots,X_{d}^{\varepsilon})$ of the follicle number in each maturity stage, the functions $\lambda_{i},\mu_{i}$, for $i\in\{0,\cdots,d\}$ are given by Eq. (4) and (9) and $\left(P_{i}^{+},P_{i}^{-}\right)_{i\in\{0,\cdots,d\}}$ are independent standard Poisson processes. In Eq. (40), $X_{0}^{\varepsilon,in}$ is an $\varepsilon\mathbb{N}$-valued random variable, and each $X_{i}^{\varepsilon,in}$, $i\in\{1,\dots,d\}$, is a $\mathbb{N}$-valued random variable.

SDE (40) defines uniquely (in law) a continuous time Markov chain in $\varepsilon{\mathbb{N}}\times{\mathbb{N}}^{d}$. We note the following conservation law,

[TABLE]

In the following, we consider a sequence $X^{\varepsilon}$ of solutions to system (40) in the limit $\varepsilon$ tends to [math].

3.2. Limit model

Formally, setting $\varepsilon=0$ in system (40) leads to the following system for $(\bar{x}_{0},\bar{f}_{t})$, coupling the dynamics of a deterministic continuous function on $\mathbb{R}_{+}$, $\bar{x}_{0}$, with those of a time-dependent measure on $\mathbb{N}^{d}$, $\bar{f}_{t}$, for all $t\geq 0$,

[TABLE]

In system (42), $x^{in}_{0}$ is a real positive constant, and, for any $x_{0}>0$, $\bar{A}_{x_{0}}$ is an operator, defined, for all bounded functions $\psi$ on ${\mathbb{N}}^{d}$ and for all $x=(x_{1},\cdots,x_{d})\in{\mathbb{N}}^{d}$, by

[TABLE]

where, for $i\in\{1,\dots,d\}$, $e_{i}$ is a unit vector of ${\mathbb{N}}^{d}$, with coordinate $1$ in the $i^{th}$ position and zero elsewhere, and $e_{d+1}$ is the null vector.

Remark 4**.**

Although the limit system (42) may appear quite different from its deterministic counterpart (11), it has the same flavour: the fast variable is in a “quasi-equilibrium” at any time $t$. Its law $f_{t}$ thus needs to solve the equilibrium of the Kolmogorov equations associated with SDE (40), which are written here with the help of the infinitesimal generator of the fast variable to ease the notations. More precisely, let $A^{\varepsilon}$ be the infinitesimal generator associated with the process $X^{\varepsilon}$, solution of (40), then for all functions $\phi:\mathbb{R}_{+}\times\mathbb{N}^{d}\to\mathbb{R}$ bounded and independent of the first variable ( $\forall(x_{0},x)\in\mathbb{R}_{+}\times\mathbb{N}^{d},\quad\phi(x_{0},x)=\psi(x)$), we have

[TABLE]

Thus $\overline{f}_{t}$ is the stationary solution corresponding to $A^{\varepsilon}$ when the slow variable is “frozen”.

System (42) is not necessarily well-posed, as there may be several solutions $\bar{f}_{t}$ for a given $\bar{x}_{0}$. In the next two specific examples, we can prove that system (42) does admit a single solution, which is a natural limit candidate for the sequence $X^{\varepsilon}$.

Example 5** (Linear case).**

Let us suppose that $K_{1,i}=K_{2,i}=0$ for all $i\in\{0,\dots,d\}$, and $(f_{i}+g_{i})>0$ for all $i\in\{1,\dots,d\}$. Then, system (42) becomes linear and has a unique solution. The invariant measure $\bar{f}_{t}$ has a product measure form

[TABLE]

with $\bar{f}^{i}_{t}$ a Poisson law on $\mathbb{N}$ of mean parameter $p_{i}x_{0}(t)$, with

[TABLE]

System (42) then reduces to

[TABLE]

It is classical that stationary distributions associated with the generator (43) are of product form for constant coefficients $\lambda_{i},\mu_{i}$ (see for instance [9, 12, 11]). Taking the product form in Eq. (44) for granted, the following calculus shows that each marginal distribution has to be a Poisson law.

Indeed, for a function $\psi$ which depends on the first variable only

[TABLE]

and for $\bar{x}_{0}>0$, we obtain, using expression (43),

[TABLE]

Hence, the solution $\bar{f}_{t}$ is such that for any bounded function $\psi_{1}$

[TABLE]

This holds in particular for $\psi_{1}(x)=\textbf{1}_{x=n}$ for any $n\in\mathbb{N}$, so that we obtain after calculus

[TABLE]

and then, by same arguments

[TABLE]

where $p_{i}$ is defined in Eq. (45).

Example 6** (A single feedback onto the quiescent follicle death rate).**

Let us suppose that $K_{1,i}=0$ for all $i\in\{0,\dots,d\}$ and $K_{2,i}=0$ for $i\in\{1,\dots,d\}$ but $K_{2,0}>0$. Then, system (42) can be simplified as

[TABLE]

where $p_{i}$ is defined in Eq. (45), and has a unique solution. The justification of system (47) follows that of Example (5) for the measure $f_{t}$, which is not directly modified by the feedback term onto the quiescent follicle death rates.

3.3. Convergence in the linear case

In this paragraph, we assume that $K_{1,i}=K_{2,i}=0$ for all $i\in\{0,\dots,d\}$ as in Example 5, and we assume further that the initial condition $X^{\varepsilon,in}=\left(X_{0}^{\varepsilon,in},X_{1}^{\varepsilon,in},\cdots,,X_{d}^{\varepsilon,in}\right)$ of system (40) is such that $\frac{1}{\varepsilon}X_{0}^{\varepsilon,in}$, $X_{1}^{\varepsilon,in}$,…, $X_{d}^{\varepsilon,in}$ are independent Poisson random variables of mean respectively $x_{0}^{in}/\varepsilon$, $x_{1}^{in}$,…, $x_{d}^{in}$. Denoting by $G^{in}$ the probability generating function of the integer-valued random vector $\left(\frac{1}{\varepsilon}X_{0}^{\varepsilon,in},X_{1}^{\varepsilon,in},\cdots,,X_{d}^{\varepsilon,in}\right)$, we thus have

[TABLE]

In such a case, one can solve explicitly system (40) for each $\varepsilon>0$. We briefly sketch the formal arguments and computations (see [9] for details). We define the probability generating function of the vector $\left(X_{0}^{\varepsilon}(t)/\varepsilon,X_{1}^{\varepsilon}(t),\cdots,,X_{d}^{\varepsilon}(t)\right)$, for $\textbf{z}\in\mathbb{R}^{d+1}$, by

[TABLE]

The infinitesimal generator $\tilde{A}^{\varepsilon}$ of the process $\tilde{X}^{\varepsilon}(t):=\left(X_{0}^{\varepsilon}(t)/\varepsilon,X_{1}^{\varepsilon}(t),\cdots,,X_{d}^{\varepsilon}(t)\right)$, is given, for all bounded functions $\psi$ on ${\mathbb{N}}^{d+1}$ and for all $n=(n_{0},n_{1},\cdots,n_{d})\in{\mathbb{N}}^{d+1}$, by

[TABLE]

Using the Kolmogorov backward equation,

[TABLE]

with $\psi(n)=\prod_{j=0}^{d}z_{j}^{n_{j}}$ (a truncation procedure is required to deal rigorously with such test functions), we obtain, using linearity of expectation, a first-order partial differential equation on $G^{\varepsilon}$ , given by

[TABLE]

It turns out that the unique solution to (52) is given by

[TABLE]

where $(x_{i}^{\varepsilon}(t))_{i=0,\cdots,d}$ is solution of the very same ODE as the linear version of Eq. (8) (whose solutions are given by Eq. (14)), with initial condition $x_{i}^{\varepsilon}(t=0)=x_{i}^{in}$, for all $i=0...d$. Thus, at any time $t$, $X_{0}^{\varepsilon}(t)/\varepsilon$, $X_{1}^{\varepsilon}(t)$,…, $X_{d}^{\varepsilon}(t)$ are independent Poisson random variables of mean (respectively) $\frac{x_{0}^{\varepsilon}(t)}{\varepsilon}$, $x_{1}^{\varepsilon}(t)$, …, $x_{d}^{\varepsilon}(t)$.

Proposition 3**.**

Assume that $K_{1,i}=K_{2,i}=0$ for all $i\in\{0,\dots,d\}$, that $f_{i}+g_{i}$ is strictly positive for $i\in\{1,\dots,d\}$ and that the initial condition $X^{\varepsilon,in}$ is such that Eq. (48) holds. Then, for all $t>0$, $X_{0}^{\varepsilon}(t)$ converges (in law) towards the deterministic value $\bar{x}_{0}(t)=x_{0}^{in}\exp(-(f_{0}+g_{0})t)$, and, for all $i=1,\cdots,d$, $X_{i}^{\varepsilon}(t)$ converges (in law) towards a Poisson random variable of mean $\bar{x}_{i}(t)=p_{i}\bar{x}_{0}(t)$, where $p_{i}$ is given by Eq. (45).

Proof.

The proof is a direct consequence of the explicit solution (53) for the probability generating function $G^{\varepsilon}$ combined with Proposition 1. As $X_{0}^{\varepsilon}(t)/\varepsilon$ is a Poisson random variable of mean $\frac{\bar{x}_{0}(t)}{\varepsilon}$, it is clear that $\mathbb{E}\left[X_{0}^{\varepsilon}(t)\right]=\bar{x}_{0}(t)$ and $var\left(X_{0}^{\varepsilon}(t)\right)=\varepsilon\bar{x}_{0}^{2}(t)\leq\varepsilon x_{0}^{in}\to 0$ as $\varepsilon\to 0$, which implies that $X_{0}^{\varepsilon}(t)$ converges in law towards $\bar{x}_{0}(t)$. For all $i=1,\cdots,d$, and $t>0$, as $x_{i}^{\varepsilon}(t)$ converges to $p_{i}\bar{x}_{0}(t)$ as \varepsilon\to{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}0}, $X_{i}^{\varepsilon}(t)$ converges (in law) towards a Poisson random variable of mean $p_{i}\bar{x}_{0}(t)$. ∎

Remark 5**.**

The slow variable $X_{0}^{\varepsilon}$ is independent of the fast variables and its dynamics are reduced to a “death and death” process with constant rate. Then convergence of this variable is given by Theorem 8.1 of [13], which enables us to obtain stronger results with fewer hypotheses on the initial condition. Suppose that

[TABLE]

then

[TABLE]

3.4. Numerical convergence

In this paragraph, we illustrate the convergence of $(X_{0}^{\varepsilon},X_{1}^{\varepsilon},\cdots,X_{d}^{\varepsilon})$ as $\varepsilon\to 0$. The chosen scenario and the parameter values are detailed in the Appendix (section 4).

In Figure 6, we plot the empirical mean trajectories (computed over $10^{4}$ sampled trajectories) in each maturity compartment ($d=10$) for the rescaled model on a time horizon $t\in(0,1)$ in the nonlinear scenario (K_{1{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},}0}>0), together with the trajectories of the analogous ODE rescaled system (8) and its limit (11). We observe that, for each compartment, the empirical mean of the SDE seems to converge to a limit value, yet this limit does not superimpose with the ODE limit solution (which is expected in a nonlinear scenario, as the ODE and SDE limits are different).

In Figure 7, using the same parameters as in Figure 6, we display the maturity distribution in the growing follicle population, for various $\varepsilon$. We empirically quantify the convergence rate using the following error, at time $t$,

[TABLE]

and the cumulative error on time interval $[0,T]$,

[TABLE]

which can be assessed numerically as

[TABLE]

where $t_{k}=k\delta_{t}$, for k=0{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\cdots}N_{T}. In practice, we also replace the limit model $\bar{X}$ by the numerically evaluated limit model $X^{\varepsilon}$ with $\varepsilon=0.001$. We then observe that the error decreases roughly linearly with $\varepsilon$.

In Figure 8, we use the linear scenario (K_{1{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},}0}=0) detailed in the Appendix (section 4) to visualize the convergence of the fast variable of the rescaled SDE to the “quasi-stationary” distribution of the limit model. The marginals of the rescaled model are evaluated over $10^{4}$ sampled trajectories at time $t=0.1$ and $t=1$. The errors between the marginal laws of the rescaled and limit models are quantified by the total variation (restricted on the support of the numerically assessed limit model):

[TABLE]

where

[TABLE]

and $\pi_{X}$, $\pi_{y}$ are, respectively, the laws of $X$ and $Y$. The error seems to decrease in a sub-linear manner with $\varepsilon$, with a plateau for $\varepsilon<10^{-2}$, which is probably due to the limited finite sampling size ($10^{4}$).

4. Appendix - parameter values

In the numerical illustrations provided throughout the previous sections, we refer to either a linear or nonlinear scenario. As far as parameter values, the only difference is that parameter $K_{1,0}$ is set to [math] in the linear scenario. All other parameters are identical and chosen as explained below.

Numerical simulation of the PDE model

We begin by shaping the desired solution $H(x)$ and we choose functions $f$ and $g$ accordingly:

[TABLE]

[TABLE]

Motivated by our biological application, we more specifically select a two-bump function

[TABLE]

with $s=0.1$, $x_{1}=0.2$, $p_{1}=0.7$, $x_{2}=0.7$ and $p_{2}=0.3$.

Except $K_{1,0}$, all coefficients weighting the nonlinear terms, $K_{20}$, $K_{1}$, $K_{2}$, are set to zero. As a result, functions $b$, $\omega_{1}$ and $\omega_{2}$, that were introduced for the sake of genericity, are not used in the numerical illustrations.

The basal activation rate $f_{0}$ is set to $f_{0}=1$, and the basal death rate in quiescent follicles $g_{0}$ is set to $g_{0}=0.1$.

Since the population feedback onto the activation rate is mainly exerted by follicles in an intermediate maturity stage, we choose $a(x)=1_{[0.3,0.7]}$, given that the state space lies in $x\in[0,1]$. The feedback gain is set to $K_{1,0}=2$.

Finally, the time horizon covers $t\in(0,1)$, and the initial condition is given by $\rho_{0}^{ini}=100$ and $\rho^{ini}\equiv 0$.

The parameter values are summed up in Table 1 and illustrated on Figure 9.

Numerical simulation of the ODE/SDE model

For the ODE- and SDE-based models, we first set the number of compartments, $d=10$, and define $x_{i}=i/d$. Then, we adapt the functions selected in the continuous PDE case and set, for $i\in\{1,d\}$,

• •

$a_{i}=a(x_{i})$

• •

$f_{i}=f(x_{i})/d$

• •

$g_{i}=g(x_{i})$

• •

$K_{1,i}=K_{2,i}=0$

All other parameters are kept as in Table 1, while the initial condition is chosen as $Y_{0}^{ini}=100$ and $Y_{i}^{ini}=0$, $i\in\{1,d\}$.

To simulate the ODE model, we use the standard python scipy.odeint, while, to simulate the SDE model, we use an exact stochastic simulation algorithm (Gillespie).