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

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.


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 (λ i ) and exit (µ 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.
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 (ε > 0) and (ii) the model in the limit ε → 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.

Initial model
Starting from the schematic model (see Eq. (1)), we formalize a system of nonlinear ordinary differential equations (ODE) as follows. Let d ∈ N * and y = (y 0 , . . . , y d ) be a function such that, for all i ∈ {0, . . . , d}, y i : t ∈ R + → 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λ 0 andμ 0 , respectively. Follicles in the intermediate compartments (1 ≤ i ≤ d − 1) are the growing follicles, and may either mature and go to the next maturation stage i + 1, at rate λ i , or die at rate µ i . Follicles in the last compartment (i = d) are named the mature follicles and can only die at rate µ d , i.e. λ d = 0 (death in this compartment corresponds to either degeneration or ovulation). All the rates (µ i , λ i ) may depend on the growing and mature follicles population (non-local interactions), which leads to the following nonlinear ODE system dy 1 (t) dt =λ 0 (y(t))y 0 (t) − λ 1 (y(t)) + µ 1 (y(t)) y 1 (t), where, for i = 0,λ 0 (y) =f with non-negative parameter constantsf 0 ,ḡ 0 , with non-negative parameter constants, 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 ∈ 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, which shows that any follicle sub-population y i is globally bounded.

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 ε 1, such that with non-negative constants f 0 , g 0 and positive initial data x in 0 , independent of ε. Note that the initial flow f 0 y in 0 = f 0 x in 0 is preserved. We then define the rescaled solution x = (x 0 , . . . , x d ) by, for all t ≥ 0, x 0 (t) = εy 0 (t/ε), and for all i ≥ 1 , x i (t) = y i (t/ε).
Then x is solution of the following system with initial condition given by x 0 (t = 0) = x in 0 , and x i (t = 0) = x in i := y in i , for i ∈ {1, . . . , d}, and where, for and λ i and µ i , for i ∈ {1, . . . , d}, are defined in Eq. (4). We note that the conservation law (5) becomes We now consider the limit for which the small parameter ε tends to 0 and the associated sequence (x ε ) solution of system (8). In such a case, system (8) is called a "slow-fast" system (x ε 0 is the slow variable, (x ε 1 , . . . , x ε d ) are the fast variables) and the study of the limit behavior when ε → 0 is a singular perturbation problem (see for instance [21]).

Limit model
Formally, setting ε = 0 in system (8) leads to the following system: with initial condition given byx 0 (t = 0) = x in 0 , andx i (t = 0) =x in i ≥ 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 (x 1 , . . . ,x d ) can be seen as reaching instantaneously (at any time t) a quasi-steady state, "driven" by the time-dependent variablex 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 ε .

Convergence in the linear case
In this paragraph, we assume that K 1,i = K 2,i = 0 for all i ∈ {0, . . . , d}, and f i + g i > 0 for all i ∈ {1, . . . , d} as in Example 1. In such a case, one can solve system (8) explicitly for each ε > 0. The solution is given by, using vectorial notations, for all t ≥ 0, where e 1 = (1, 0, · · · , 0) T and It is thus clear that x ε 0 =x 0 is a constant sequence in ε (as both x ε 0 andx 0 have same initial conditions, and same evolution equation). For the fast variables, we prove the following where (x 1 , · · · ,x d ) is given in Eq. (12).
Proof. From Eq. (14) and initial condition, and using integration by parts, we obtain As x ε 0 is uniformly bounded (both in ε and time) by x in 0 , we obtain, taking the 1-norm, The third inequality above was deduced from Follicule number Follicule number We verify that for all t ≥ 0, Then, we obtain, which converges to 0 as ε converges to 0.
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 0 (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 , · · · , x ε d )(t) is not uniformly bounded as ε → 0.

Numerical convergence
In this paragraph, we illustrate the convergence of (x ε 0 , x ε 1 , · · · , x ε d ) towards (x 0 ,x 1 , · · · ,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 ∈ (0, 1). In each compartment, the trajectories of the rescaled model get closer and closer to the limit model as ε → 0. For ε = 0.001, they are almost indistinguishable. Note however that, for

Error
Convergence rate , for the rescaled variables x ε i , for different ε (solid colored lines, see legend insert) and the reduced limit variablesx i (black dashed lines). On the right panel, we plot the discrete l 1 norm error E 1 (t) at a fixed time t = 0.1 and t = 1 (solid red and green lines, resp.) and the l 1 -cumulative error E 2 (1) on t ∈ (0, 1) (black solid line) as a function of ε (see details in the text). The black dashed line is the straight line of slope 1 according to ε. See the Appendix (section 4) for details on the parameter values used in the numerical simulations.
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 ε. 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, and the cumulative error on time interval [0, T ], which can be assessed numerically asẼ where t k = kδ t , for k = 0· · ·N T . From the right panel of Figure 2, we can see that, at least for a small enough ε, the error is inversely proportional to ε, E i ≈ C te ε −1 , where the constant pre-factor may depend on the chosen norm or the particular time t.

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.

Rescaled model
Denoting by ρ 0 (t) the number of quiescent follicles and by ρ(t, x) the population density of follicles of maturity x, we consider the following coupled ODE-PDE system, for all t ≥ 0, where and, for all x ∈ (0, 1), with initial condition We assume that f, g, K 1 , K 2 , a, b, w 1 , w 2 , ρ 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,5,6] and references therein), yet this is beyond the scope of this work. We can write the following conservation law, that gives (at least formally) In the following, we consider a sequence (ρ ε 0 , ρ ε ) of solutions of system (20) in the limit ε → 0.

Limit model
Formally, setting ε = 0 in system (20) leads to the following system: for all t ≥ 0, with an initial condition given byρ 0 (t = 0) = ρ in 0 , andρ(t = 0, .) =ρ 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ρ for a givenρ 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 (ρ ε 0 , ρ ε ).

Convergence in the linear case
In this paragraph, we assume that K 1,0 = K 2,0 = 0, K 1 ≡ 0, K 2 ≡ 0, and f (0) > 0 as in Example 3. In such a case, one can solve explicitly system (20) for each ε > 0 using the characteristics method. We obtain where X(0; t, x) is the location of the characteristic at time 0, given that it goes through the point x at time t, namely: It is thus clear that ρ ε 0 =ρ 0 is a constant sequence in ε. For the population density ρ ε , which is here the fast unknown, we prove the following whereρ is given in Eq. (26).
Proof. It is clear that for any η > 0, there exists ε such that for all ε < ε , and all t > η we have Then, comparing the solutions of Eq. (29) and Eq. (26) allows us to conclude that, for all t > η, Remark 3. As in Remark 1, we can see that during a time of order εt, we cannot get the convergence of the rescaled model towards the reduced one, which precludes uniform convergence in time starting from t =0.

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.

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 λ 0 , λ and µ upon the solutionρ(t, x), which itself depends onρ 0 . We propose to treat this nonlinearity with a fixed point scheme. At each time step t n = n∆ t , for n = 0, . . . , N with T = N ∆ t we build a sequence ρ (t n , x) such that lim →∞ ρ (t n , x) =ρ(t n , x).

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ρ 0 (t), and we use Eq. (27) to computeρ(t, ·). We use this reference solution to assess the performances of the numerical scheme described in paragraph 2.4.1. We display on In the top right panel, we display the relative errors in L 1 and L ∞ norms between the pseudo exact and the numerical solutionρ 0 (t), as a function of ∆ t = ∆ x = 1/N : ..,N |ρ 0 (t n ) −ρ n | max n=0,...,N |ρ 0 (t n )| , and the relative errors between the pseudo exact and the numerical solutionsρ(T, x)

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∆ t , for n = 0, . . . , and, for x k = k∆ x , k = 0, . . . , M , with M ∆ x = 1, We integrate numerically the PDE between t n and t n+1 and over [x k , x k+1 ] by freezing the nonlinear coefficients λ(ρ(t, .), x) and µ(ρ(t, .)) at time t n ρ ε,n+1 At each time step, we compute both the PDE and ODE coefficients using the midpoint rule and the numerical solution (ρ ε,n k ) k=0,...,M as a piecewise constant solution For the explicit scheme (32), two stability conditions must be satisfied • CFL-like stability condition: which can be rewritten as (34) • Positivity conservation condition: 1+ which can be rewritten as ∆ n t ε∆ x ≤ 1 max k (λ ε,n k + ∆ x µ ε,n k ) .
The overall numerical scheme proceeds as follows: (1) Initialization : ρ ε,0 k = ρ in (x k ) andρ ε,0 0 = ρ in 0 .  Figure 4. In the left panel, we see that ρ 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.

ε-convergence towards the limit model
So far we have proved the convergence of the rescaled model towards the limit model when ε → 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 ε-convergence of (ρ ε 0 , ρ ε ) towards (ρ 0 ,ρ) 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 (ρ 0 ,ρ) are simulated using Eq. (28) and Eq. (27). The results are displayed in Figure 5 for a set of ε values and two discretizations M = N = 100 and M = N = 200. The agreement of ρ 0 (t) with the limit model solution (top left panel), and that of ρ ε (t, x) (bottom left panel), are qualitatively good as soon as ε ≤ 10 −2 . The relative error for ρ 0 (t) (top right panel) exhibits a linear behavior in ε. The error curves for the convergence of the solution ρ ε (t, x) in the domain (bottom right panel) are not linear, and remain beyond a threshold when ε goes to 0. However, both the value of ε 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 ε. Since our current numerical scheme is explicit in time, any refinement must be simultaneous in ∆ x and ∆ t , and as a consequence the cost in CPU time depends quadratically in ε −1 . In practice, realistic values for ε 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.

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.

Rescaled model
We consider the following coupled Poisson-driven SDE system, given by, for all t ≥ 0, where X ε represents the vectorial process (X ε 0 , · · · , X ε d ) of the follicle number in each maturity stage, the functions λ i , µ i , for i ∈ {0, · · · , d} are given by Eq. (4) and (9) and P + i , P − i i∈{0,··· ,d} are independent standard Poisson processes. In Eq. (36), X ε,in 0 is an εN-valued random variable, and each X ε,in i , i ∈ {1, . . . , d}, is a Nvalued random variable. SDE (36) defines uniquely (in law) a continuous time Markov chain in εN × N d . We note the following conservation law, In the following, we consider a sequence X ε of solutions to system (36) in the limit ε tends to 0.

Limit model
Formally, setting ε = 0 in system (36) leads to the following system for (x 0 ,f t ), coupling the dynamics of a deterministic continuous function on R + ,x 0 , with those of a time-dependent measure on N d ,f t , for all t ≥ 0, In system (38), x in 0 is a real positive constant, and, for any x 0 > 0,Ā x0 is an operator, defined, for all bounded functions ψ on N d and for all x = (x 1 , · · · , x d ) ∈ N d , bȳ where, for i ∈ {1, . . . , d}, e i is a unit vector of 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 (38) 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 (36), which are written here with the help of the infinitesimal generator of the fast variable to ease the notations. More precisely, let A ε be the infinitesimal generator associated with the process X ε , solution of (36), then for all functions φ : R + × N d → R bounded and independent of the first variable ( Thus f t is the stationary solution corresponding to A ε when the slow variable is "frozen". System (38) is not necessarily well-posed, as there may be several solutionsf t for a givenx 0 . In the next two specific examples, we can prove that system (38) does admit a single solution, which is a natural limit candidate for the sequence X ε .
Example 5 (Linear case). Let us suppose that K 1,i = K 2,i = 0 for all i ∈ {0, . . . , d}, and (f i + g i ) > 0 for all i ∈ {1, . . . , d}. Then, system (38) becomes linear and has a unique solution. The invariant measuref t has a product measure form withf i t a Poisson law on N of mean parameter p i x 0 (t), with System (38) then reduces to It is classical that stationary distributions associated with the generator (39) are of product form for constant coefficients λ i , µ i (see for instance [9,11,12]). Taking the product form in Eq. (40) for granted, the following calculus shows that each marginal distribution has to be a Poisson law. Indeed, for a function ψ which depends on the first variable only and forx 0 > 0, we obtain, using expression (39), Hence, the solutionf t is such that for any bounded function ψ 1 Example 6 (A single feedback onto the quiescent follicle death rate). Let us suppose that K 1,i = 0 for all i ∈ {0, . . . , d} and K 2,i = 0 for i ∈ {1, . . . , d} but K 2,0 > 0. Then, system (38) can be simplified as where p i is defined in Eq. (41), and has a unique solution. The justification of system (43) 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.
Remark 5. The slow variable X ε 0 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

Numerical convergence
In this paragraph, we illustrate the convergence of (X ε 0 , X ε 1 , · · · , X ε d ) as ε → 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 ∈ (0, 1) in the nonlinear scenario (K 1,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 ε. We empirically quantify the convergence rate using the following error, at time Follicule number .001 ODE Limit model Figure 6. In the same way as in Figure 1, we plot the trajectories in each maturity compartment (d = 10), for the rescaled variables X ε i of the SDE system (36) (solid lines) and the ODE system (8) (dashed lines), for different ε (see legend). For the SDE, we plot the empirical mean computed over 10000 trajectories. The limit variableX i of the ODE system corresponds to the black dashed line. t, and the cumulative error on time interval [0, T ], which can be assessed numerically as where t k = kδ t , for k = 0· · ·N T . In practice, we also replace the limit modelX by the numerically evaluated limit model X ε with ε = 0.001. We then observe that the error decreases roughly linearly with ε.
In Figure 8, we use the linear scenario (K 1,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

Error
Convergence rate O(ε) ε = 10 0.00 ε = 10 −0.33 ε = 10 −0.67 ε = 10 −1.00 ε = 10 −1.33 ε = 10 −1.67 ε = 10 −2.00 ε = 10 −2.33 ε = 10 −2.67 ε = 10 −3.00 ODE limit model , for the rescaled variables X ε i , for different ε (solid colored lines, see legend insert) and the limit variableX i (black dashed line). On the right panel, we plot the discrete l 1 norm error E 1 (t) at a fixed time t = 0.1 and t = 1 (solid red and green lines, resp.) and the l 1cumulative error E 2 (1) on t ∈ (0, 1) (black solid line) as a function of ε (see details in body text). The black dashed line is the straight line of slope 1 according to ε. 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): where d T V (X, Y ) = max{| π X (i) − π Y (i) | , i ∈ N} , and π X , π y are, respectively, the laws of X and Y . The error seems to decrease in a sub-linear manner with ε, with a plateau for ε < 10 −2 , which is probably due to the limited finite sampling size (10 4 ).  . Empirical law of X ε i in each maturity compartment at time t = 1 for different ε in colored bars (see legend insert) and the limit distributionX i (black solid lines). On the bottom leftmost panel, we plot the total variation error E(t) at a fixed time t = 0.1 and t = 1 (solid red and green lines, resp.). The black dashed line is the straight line of slope 1 according to ε.

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 0 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: 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, ω 1 and ω 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 ∈ [0, 1]. The feedback gain is set to K 1,0 = 2. Finally, the time horizon covers t ∈ (0, 1), and the initial condition is given by ρ ini 0 = 100 and ρ ini ≡ 0. The parameter values are summed up in Table 1 and illustrated on Figure 9.
• K 1,i = K 2,i = 0 All other parameters are kept as in Table 1, while the initial condition is chosen as Y ini 0 = 100 and Y ini i = 0, i ∈ {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).