# Bifurcation analysis of a size-structured population model: Application to oocyte dynamics and ovarian cycle

Frédérique Clément, Louis Fostier, Romain Yvinec

► To cite this version:

Frédérique Clément, Louis Fostier, Romain Yvinec. Bifurcation analysis of a size-structured population model: Application to oocyte dynamics and ovarian cycle. SIAM Journal on Applied Dynamical Systems, 2025, 24 (Issue 3), pp.2427 - 2472. 10.1137/24M1705147. hal-04699357v2

HAL Id: hal-04699357

https://hal.science/hal-04699357v2

Submitted on 5 Dec 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# BIFURCATION ANALYSIS OF A SIZE-STRUCTURED POPULATION MODEL: APPLICATION TO OOCYTE DYNAMICS AND OVARIAN CYCLE *

FRÉDÉRIQUE CLÉMENT†, LOUIS FOSTIER†‡, AND ROMAIN YVINEC†‡

**Abstract.** We introduce and analyze a quasilinear size-structured population model with non-linearities accounting for nonlocal interactions between individuals. The recruitment (immigration), growth and death rates are inhomogeneous in time and/or space and depend on weighted averages of the density. We first prove the existence and uniqueness of globally bounded weak solutions using the characteristic curves and Banach fixed point Theorem, after transforming the partial differential equation into an equivalent system of integral equations. We then investigate the long-time behavior of the PDE in the case when the growth rate is separable. Applying a classical time-scaling transformation, the problem boils down to a PDE with linear growth rate and nonlinear inflow boundary condition, entering the theoretical framework of abstract semilinear Cauchy problems. We can then perform a bifurcation analysis which reveals the richness of the model behavior. Depending on the ratio of the recruitment to the growth rate, the model can exhibit multistability and stable oscillatory solutions, emanating respectively through saddle-node and Hopf bifurcations. We illustrate these theoretical results on the biological application motivating this work, oogenesis, the process of production and maturation of female gametes (oocytes) that is critical to reproductive fitness.

**Key words.** size-structured population dynamics; global existence; asymptotic behavior; bifurcation analysis; reproductive biology

MSC codes. 35B32; 35B40; 35Q49; 92D25

## 1. Introduction.

### 1.1. Model and state of the art.
We consider the following quasilinear initial boundary value problem modeling the evolution of a continuous size-structured population with nonlocal interactions between individuals, over a bounded domain:

$$ (1.1) \qquad \left\{ \begin{array}{l} \partial_t \rho(t, x) + \partial_x (\Lambda(P(t), x)\rho(t, x)) = -\mu(P(t), x)\rho(t, x), \quad x \in (0, 1), \ t > 0 \\ \Lambda(P(t), 0)\rho(t, 0) = r(t, P(t)), \quad t > 0 \\ \rho(0, x) = \rho_0(x), \quad x \in (0, 1). \end{array} \right. $$

Here $\rho(t, x)$ stands for the density of individuals at time $t$ with respect to size $x$, $\Lambda$ is the individual growth rate and $\mu$ the death rate. In the boundary condition (second line of Eq. (1.1)), $r$ is classically called the birth term, but represents in our study the inflow of new individuals into the population (similar to an immigration process) and is termed the *recruitment rate*. The *interaction terms*

$$ (1.2) \qquad P(t) = (P_i(t))_{1 \le i \le N}, \quad P_i(t) = \int_0^1 \omega_i(x) \rho(t, x) dx, $$

<footnote>†Université Paris-Saclay, Inria, Centre Inria de Saclay, 91120, Palaiseau, France (fred-erique.clement@inria.fr).</footnote>

<footnote>‡PRC, INRAE, CNRS, Université de Tours, 37380 Nouzilly, France (louis.fostier@inria.fr, ro-main.yvinec@inria.fr).</footnote>
are weighted averages of the density $\rho$ at time $t$, where the weight functions $\omega := (\omega_i)_{1\le i\le N}$ select (smooth) subpopulations amounting to size classes that contribute to the interaction term $P(t)$. Generally, both the growth rate $\Lambda$ and the death rate $\mu$ depend on size $x$ and the interaction term $P(t)$, while the recruitment rate depends on time $t$ and the interaction term $P(t)$. When the growth rate $\Lambda$ is constant, Eq. (1.2) models an age-structured population. If $\Lambda$ does not depend on $P$, the problem is semilinear. It is worth mentioning that a recruitment term expressed as the sum of an inflow of new individuals $C(t, P(t))$ and a birth term with a birth rate of the form $\beta(t, P(t), x) = \tilde{\beta}(t, P(t))\omega_1(x)$ can be seen as a specific case of Eq. (1.1)-(1.2):

$$r(t, P(t)) = C(t, P(t)) + \int_{0}^{1} \beta(t, P(t), x)\rho(t, x)dx = C(t, P(t)) + \tilde{\beta}(t, P(t))P_{1}(t).$$

Structured population models, which describe population dynamics in terms of individual characteristics such as age or size, have been the matter of extensive studies (see [47, 38, 27, 13, 37] for a good overview). They have found concrete applications in diverse fields such as cell population dynamics [25, 16, 23, 6], ecology [44, 43], or manufacturing system modeling [42, 12, 33].

The well-posedness of quasilinear size-structured models has been extensively studied. The existence of solutions to these models is typically demonstrated through two approaches: the characteristic method combined with a fixed-point argument, such as the Banach Theorem [45, 9, 30], or the Schauder Theorem [31](the latter ensuring only local existence but under weaker assumptions), or by finite difference approximation [1]. The model formulated in Eq. (1.1)-(1.2) differs slightly from previous works, as it involves a positive growth rate at $x = 1$ (leading to positive outflow) and a recruitment rate $r$ not explicitly formulated as a birth term. However, the proof for the global existence-uniqueness of solutions to Eq. (1.1)-(1.2) is based on the Banach fixed-point Theorem under standard assumptions (given in Sec. 2.1).

The long-time behavior of solutions to quasilinear size-structured population models is challenging to study. While for semilinear population models the Principle of Linearized Stability for stationary solutions [47, 29, 24], and Hopf bifurcation Theorem [36, 34, 37] (associated with oscillatory solutions) have been established, this is not the case for general quasilinear size-structured models. Recent works, such as those by Diekmann et al. [8, 5, 22], focus on formalizing the relationship between the PDE formulation and the renewal equation, for which long-time behavior results like the Principle of Linearized Stability are already known [14, 15, 5].

In the particular case when a separable growth rate $\Lambda(P, x) = g(P)f(x)$ is assumed, a classical implicit time-scaling transformation can remove the nonlinearity of the growth rate [24, 25, 10], so that the long-time behavior of the quasilinear model can be studied using results from semilinear equation theory.

Regarding long-time behavior, Farkas et al. have investigated the stability of linearized PDEs around stationary solutions in various quasilinear and semilinear models that explicitly include a recruitment term composed of an external inflow of individuals and a birth term. Specifically, they analyzed a general size-structured model with constant inflow and nonlinear birth rate [20], a general age-structured model with nonlinear (density-dependent) inflow and linear birth rate [19], and a classical Daphnia model with nonlinear (food-dependent) external inflow and nonlinear birth rate [21]. Notably, the bistability result in [20] aligns with our findings in Sec. 4.2.

The current work focuses on the long-time behavior of Eq. (1.1) with separable growth
speed, no death, and a recruitment term that does not explicitly include a birth term. As a result, the existence and stability of stationary solutions can be analyzed in a more straightforward way than in a general size-structured model, and, as said before, the Principle of Linearized Stability can be applied thanks to the implicit time-scaling. Notably, unlike in Farkas’s studies, multiple interaction terms $P \in \mathbb{R}^N$ are considered rather than just the total mass $P = \int_0^1 \rho(t, x)dx$, and the implications of this difference will be discussed. Additionally, this work explores the existence of Hopf bifurcations and periodic solutions, which were not addressed in Farkas’s works, and, as far as we know, have not been addressed in the literature for quasilinear models with separable growth speed.

**1.2. Biological and biomathematical motivations.** The main biological application motivating the study of Eq. (1.1) lies in developmental/reproductive biology. We are interested in oogenesis, the process of growth and maturation undergone by germ cells until the release of a fertilizable female gamete, with a focus on fish oogenesis and spawning cycle.

Understanding fish oogenesis is useful for basic research, fish farming and ecophysiology. Indeed, due to their aquatic lifestyle, fish are particularly exposed to pollutants such as endocrine disruptors (ED), while, due to their reproductive physiology subject to chemically-induced sex reversal and disruption of the sex-ratio, they are particularly sensitive to ED. Hence fish are sentinel vertebrate species for chemical pollution enabling one to study the multiscale impact of ED from individual reproductive troubles up to population collapse [32].

Once the germ cells have reached the oocyte stage, they stop proliferating and pursue their development enclosed in somatic structures called ovarian follicles. In fish, there is a single layer of somatic cells surrounding the oocyte, so that follicle growth is almost exclusively due to the increase in oocyte diameter. Hence oocyte size is a natural structuring variable for the whole oocyte population within the ovary, and a proxy for oocyte maturity which appears as a main entry into the reference classifications of oogenesis stages [28, 41].

The formalism of structured population dynamics with nonlocal nonlinearities is well suited for representing the controls exerted onto oogenesis, which emanate, either directly (i.e. on the local ovarian scale) or indirectly (i.e. through endocrine loops involving the hypothalamus and pituitary) from the oocyte and somatic companion cells. We previously handled such a formalism in the framework of mammal oogenesis to investigate the earliest steps of follicle development (i.e. activation from a non-renewed pool of quiescent follicles) [3, 4] and the competition between terminally developing follicles for ovulation [2].

Here, we apply the generic model formulated in Eq. (1.1) to design a size-structured model of (follicle-enclosed) oogenesis in fish, where $x$ corresponds to the oocyte size. The main salient features are (i) the non-zero recruitment rate (boundary condition at $x=0$), since, contrary to mammals, oocytes can be renewed from germline stem cells all along life, and (ii) the non-scalar formulation of the nonlocal control term, $P(t)$, to account for the combination of size-specific hormonal signals as well as for the size-dependent effect of those signals on the growth speed. The outflow (at $x=1$) corresponds to the spawning event, so that the integral in time of the outflow assesses the cumulative number of spawned eggs.

In the field of mathematical biology, our size-structured approach is actually a more generic physiologically-structured approach that could be adapted to other contexts, in which the structuring variable would monotonously increase along progres-
sion in a (cell) lineage, and the control of cell commitment and maturation could be represented as nonlocal interactions mediated by the (possibly weighted) population size. This could for instance be the case for the erythroid lineage, characterized by an increasing hemoglobin content from immature precursors to red cells [25]. The marrow-kidney feedback loop involving the synthesis of erythropoietin (EPO) from oxygen-sensing kidney cells, and then stimulation of precursors by EPO, could also be designed in a compact way, from the number of red cells, possibly weighted by their hemoglobin content.

**1.3. Outline of this work.** In Sec. 2, under some regularity assumptions, and a crucial bounding assumption on the recruitment rate, we derive a global existence-uniqueness result on the solutions to Eq. (1.1)-(1.2), and a first result on the existence of stationary solutions. In Sec. 3, with a separable growth rate and in the absence of death, we exhibit a simpler characterization of the stationary solutions, and investigate with care the influence of the recruitment-to-growth ratio on their existence as well as local stability. Using these results, we fully describe in Sec. 4 the long-time behavior of the model when the recruitment-to-growth ratio is monotonous in all interaction terms, with same direction of variation in all these terms. We first illustrate a scenario where negative feedback on the recruitment rate and positive feedback on the growth rate can lead to monostability or oscillations, with the latter occurring through a Hopf bifurcation. We then exhibit a situation where negative feedback affects both the recruitment rate and the growth rate, leading to either mono- or multistability through two successive saddle-node bifurcations. These cases are also motivated by the study of fish oocyte populations over the entire life span. In the two considered scenarios (Hopf and saddle-node bifurcations), we include numerical simulations involving a slow time-varying recruitment rate, which leads to transitions in the reproductive status of aged individuals.

**2. Global Solution and Stationary Solution Analysis.**

**2.1. Notations and assumptions.** We introduce the notations to be used consistently in this study. Let $I$ be an interval.

• The essential supremum is denoted by $\|\phi\|_{\infty,I} := \ess\sup_{y \in I} |\phi(y)|$, with $\phi \in L^{\infty}(I)$. When $I$ is obvious, we will just note $\|\phi\|_{\infty} = \|\phi\|_{\infty,I}$. In particular, $\|\omega_i\|_{\infty} = \|\omega_i\|_{\infty,[0,1]}$ for $1 \le i \le N$, and $\|\rho_0\|_{\infty} = \|\rho_0\|_{\infty,[0,1]}$.

• For any $\omega = (\omega_1, \dots, \omega_N) \in (L^\infty(I))$, we define $\|\omega\|_\infty := \sum_{i=1}^N \|\omega_i\|_\infty$.

• For any $\phi \in L^1(I)$ and $\omega \in (L^\infty(I))^N$, we note $\int_I \omega \phi := (\int_I \omega_1 \phi, \dots, \int_I \omega_i \phi, \dots, \int_I \omega_N \phi) \in \mathbb{R}^N$.

• A vector $P \in \mathbb{R}^N$ is said nonnegative (denoted by $P \ge 0$), if $P_i \ge 0$ for all $i = 1, \dots, N$, and positive (denoted by $P > 0$) if $P_i > 0$ for all $i = 1, \dots, N$. We make the following assumptions throughout the paper, unless specified otherwise:

(H1) $\rho_0 \in L^\infty(0, 1)$ is nonnegative.

(H2) $\Lambda : [0, +\infty)^N \times [0, 1] \mapsto \mathbb{R}_+^*$ is strictly positive, continuously differentiable with respect to $x$, and Lipschitz on every compact set $K$ over $[0, +\infty)^N \times [0, 1]$, with Lipschitz constant $L_{\Lambda,K}$.

(H3) $\mu : [0, +\infty)^N \times [0, 1] \mapsto \mathbb{R}_+$ is nonnegative and Lipschitz continuous on every compact set $K$ over $[0, +\infty)^N \times [0, 1]$, with Lipschitz constant $L_{\mu,K}$.

(H4) $r : [0, +\infty) \times [0, +\infty)^N \mapsto \mathbb{R}_+$ is nonnegative and Lipschitz continuous on every compact set $K$ over $[0, +\infty) \times [0, +\infty)^N$, with Lipschitz constant $L_{r,K}$.
(H5) There exist $(r_i)_{0 \le i \le N}$ locally bounded functions, i.e. in $L^\infty(0, T)$ for all $T > 0$, such that for all $t \ge 0$ and $P \ge 0$,

$$ (2.1) \qquad |r(t, P)| \le r_0(t) + \sum_{i=1}^{N} r_i(t) P_i. $$

(H6) For all $i \in [1; N]$, $\omega_i : [0, 1] \to \mathbb{R}$ is nonnegative and Lipschitz continuous, with Lipschitz constant $L_{\omega_i}$, and $L_{\omega} := \sum_{i=1}^{N} L_{\omega_i}$.

In assumption (H1), the initial distribution $\rho_0$ is considered in $L^\infty$ rather than $L^1$. Given assumptions (H2)-(H6), we ensure that the distribution remains $L^\infty$ at all times, and our stationary solution study is restricted to $L^\infty$ functions as well. Assumption (H1) is required for the contraction strategy employed to establish the existence-uniqueness Theorem 2.2. Assumption (H2) implies that the growth rate has a strictly positive infimum over any compact set, preventing mass concentration phenomena and thereby justifying to work on $L^\infty$. Inequality Eq. (2.1) in assumption (H5) ensures that not too many individuals enter the system simultaneously, which is essential for defining global-in-time solutions. By integrating Eq. (1.1) over the spatial domain, it becomes evident that a quadratic recruitment term in P can lead to a finite-time explosion of the total mass.

We do not claim that assumptions (H2)-(H6) are optimal, but they are rather generic and flexible enough to cope with a range of applications.

## 2.2. Existence and uniqueness.
From assumption (H2), we can define the classical notion of characteristic curves and the associated diffeomorphisms (see Figure 8 in Appendix A to aid in understanding the definitions and notations).

**DEFINITION - PROPOSITION 2.1.** For a given time $T > 0$, and $P$ a nonnegative continuous function from $[0, T]$ to $\mathbb{R}^N$, we define $X_P(s; t, x)$ as the characteristic curve at time $s$ passing through $x$ at time $t$. More precisely, for any $(t, x) \in [0, T] \times (0, 1)$, $X_P(s; t, x)$ is the unique maximal solution to the following ordinary differential equation on $[0, T] \times (0, 1)$:

$$ (2.2) \qquad \dot{X}(s) = \Lambda(P(s), X(s)), \quad X(t) = x. $$

Let $I_P(t, x) = \{s \in [0, T], X_P(s; t, x) \in (0, 1)\}$ be the maximal interval definition of $X_P(\cdot; t, x)$, and

$$ \alpha_P(t, x) = \inf I_P(t, x), \quad \beta_P(t, x) = \lim_{s \to \alpha_P(t, x)} X_P(s; t, x) = X_P(\alpha_P(t, x); t, x). $$

Finally, we define

$$ X_P(s; t, 1) = \lim_{x \to 1} X_P(s; t, x), \quad X_P(s; t, 0) = \lim_{x \to 0} X_P(s; t, x), \quad X_{0,P} := X_P(s; 0, 0). $$

For each $t \in [0, T]$, $\alpha_P(t, \cdot)$ defines a $C^1$-diffeormorphism from $D_\alpha(t) := (0, X_{0,P}(t))$ to $(0, t)$ and $\beta_P(t, \cdot)$ defines a $C^1$-diffeormorphism from $D_\beta(t) := (X_{0,P}(t), 1)$ to $(0, \beta_P(t, 1))$. Furthermore, on domains $D_\alpha$ (respectively $D_\beta$), the derivatives of $\alpha_P$
(respectively $\beta_P$) are given by

$$ (2.3) \qquad \begin{aligned} \frac{\partial \alpha_P}{\partial t}(t, x) &= \frac{\Lambda(P(t), x)}{\Lambda(P(\alpha_P(t, x)), 0)} e^{-\int_{\alpha_P(t, x)}^t \Lambda_x(P(\theta), X_P(\theta; t, x)) d\theta}, \\ \frac{\partial \alpha_P}{\partial x}(t, x) &= -\frac{1}{\Lambda(P(\alpha_P(t, x)), 0)} e^{-\int_{\alpha_P(t, x)}^t \Lambda_x(P(\theta), X_P(\theta; t, x)) d\theta}, \\ \frac{\partial \beta_P}{\partial t}(t, x) &= -\Lambda(P(t), x) e^{-\int_0^t \Lambda_x((P(\theta), X_P(\theta; t, x)) d\theta}, \\ \frac{\partial \beta_P}{\partial x}(t, x) &= e^{-\int_0^t \Lambda_x(P(\theta), X_P(\theta; t, x)) d\theta}. \end{aligned} $$

Given $P$, the existence and uniqueness of an explicit solution through the characteristic formula is classical (see beginning of Appendix A). The main technical challenge of this section is to show that it can be extended to the general nonlinear case when $P$ is given by Eq. (1.2).

The solutions considered here are functions that satisfy the variational formulation Eq. (2.4) of Eq. (1.1). Within this framework, one can accommodate $L^\infty$ initial conditions that may be incompatible with the boundary conditions, leading to discontinuous solutions along the first characteristic curve.

**DEFINITION 2.1.** *Considering $T > 0$, we say that $\rho$ is a weak solution to Eq. (1.1) if $\rho \in C^0([0, T]; L^1(0, 1)) \cap L^\infty((0, T) \times (0, 1))$, and for all $\tau \in [0, T]$ and $\varphi \in C^1([0, \tau] \times [0, 1])$ such that $\forall x \in [0, 1]$, $\varphi(\tau, x) = 0$, and $\forall t \in [0, \tau]$, $\varphi(t, 1) = 0$, we have:*

$$ (2.4) \qquad \begin{aligned} & \int_0^\tau \int_0^1 \rho(t, x) [\partial_t \varphi(t, x) + \Lambda(P(t), x) \partial_x \varphi(t, x) - \mu(P(t), x) \varphi(t, x)] dxdt \\ & + \int_0^\tau r(t, P(t)) \varphi(t, 0) dt + \int_0^1 \rho_0(x) \varphi(0, x) dx = 0, \end{aligned} $$

where $P$ is given by Eq. (1.2).

We now state the main result of this section, which guarantees the existence and uniqueness of a global weak solution to Eq. (1.1), satisfying a quasi-explicit characteristic formula.

**THEOREM 2.2.** Let $T > 0$. Under assumptions (H1)-(H6), Eq. (1.1) admits a unique weak solution $\rho$ in the sense of definition 2.1. Furthermore, $\rho$ verifies: (2.5)

$$ \left\{ \begin{array}{ll} \rho(t, x) = \rho_0(\beta_P(t, x)) \mathcal{F}_P(t; 0, \beta_P(t, x)) \frac{\partial \beta_P}{\partial x}(t, x), & \text{if } X_{0,P}(t) \le x \le 1, 0 \le t \le T \\ \\ \rho(t, x) = -r(\alpha_P(t, x), P(\alpha_P(t, x))) \mathcal{F}_P(t; \alpha_P(t, x), 0) \frac{\partial \alpha_P}{\partial x}(t, x) & \text{otherwise} \end{array} \right. $$

where $P$ is given by Eq. (1.2), $\alpha_P$, $\beta_P$ defined in Definition-Proposition 2.1 and $\mathcal{F}_P(s; t, x)$ is defined as the probability that an individual of size $x$ at time $t$ is still alive at time $s$, that is

$$ \mathcal{F}_P(s; t, x) := e^{-\int_t^s \mu(P(\theta), X_P(\theta; t, x)) d\theta}. $$

The principle of the proof follows Sec. 3 of [42]. Let $t > 0$ be a sufficiently small time such that $X_{0,P}(t)$ is defined. If $\rho$ is solution to Eq. (1.1), then, as shown in Proposition A.2, $\rho$ satisfies the characteristic formula Eq. (2.5) with $P$ given by Eq. (1.2).
Therefore, by direct integration, we obtain for all $i \in \{1, \dots, N\}$:

$$ 
\begin{aligned}
P_i(t) &= \int_0^1 \omega_i(y) \rho(t, y) dy \\
&= -\int_0^{X_{0,P}(t)} r(\alpha_P(t, y), P(\alpha_P(t, y))) \mathcal{F}_P(t; \alpha_P(t, y), 0) \frac{\partial \alpha_P}{\partial x}(t, y) \omega_i(y) dy \\
&\quad + \int_{X_{0,P}(t)}^1 \rho_0(\beta_P(t, y)) \mathcal{F}_P(t; 0, \beta_P(t, y)) \frac{\partial \beta_P}{\partial y}(t, y) \omega_i(y) dy.
\end{aligned}
$$

Then, by making the change of variables $u(y) = \alpha_P(t, y)$ and $x(y) = \beta_P(t, y)$ (see Eq. (2.3)), we define:

$$ (2.6) \qquad F_i(P)(t) := \int_0^t r(u, P(u)) \mathcal{F}_P(t; u, 0) \omega_i(X_P(t; u, 0)) du \\ 
\qquad \qquad + \int_0^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))) d\theta} \rho_0(x) \mathcal{F}_P(t; 0, x) \omega_i(X_P(t; 0, x)) dx $$

where $P$ must satisfy $P_i = F_i(P)$, $1 \le i \le N$.

Conversely, it can be shown (see A.4) that if $P$ is a fixed point of Eq. (2.6), then $\rho$ defined by Eq. (2.5) is indeed a weak solution to Eq. (1.1).

Thus, the key point is to prove that Eq. (2.6) admits a unique fixed point, at least locally. Our strategy is to use the Banach fixed-point Theorem for the mapping defined by Eq. (2.6). First, Lemma 2.3 allows us to control the distance between two characteristic trajectories associated with different $P$ values, leading to the contraction property established in Lemma 2.4. Next, the local unique solution is extended globally in time. The latter step requires a uniform in time control on the weighted averages of the density, which is obtained using Grönwall's Lemma and relies on the linear control on $r$ provided by assumption (H5).

We use the Banach fixed-point Theorem on the space

$$ \Omega_{\delta,M} := \{ P \in C^0([0, \delta])^N \mid \forall i \in [1, N], P_i \ge 0 \text{ and } \|P_i\|_{C^0([0, \delta])} = \sup_{t \in [0, \delta]} P_i(t) \le M \} $$

with the associated norm $\|P\|_{\Omega_{\delta,M}} := \sum_{i=1}^N \|P_i\|_{C^0([0, \delta])}$, where $\delta$ and $M$ will be chosen later. The space $\Omega_{\delta,M}$ is a product of closed sets of $C^0([0, \delta])$ which is a complete metric space. Therefore, $\Omega_{\delta,M}$ is also complete. We recall that, by (H2), $\Lambda$ is Lipschitz on $[0, M]^N \times [0, 1]$ and denote $L_{\Lambda,M} := L_{\Lambda, [0, M]^N \times [0, 1]}$ the associated Lipschitz constant.

**LEMMA 2.3.** Let $0 < \delta \le (2L_{\Lambda,M})^{-1}$. For all $P, Q \in \Omega_{\delta,M}$, $t \in [0, \delta]$, and $x \in [0, 1]$, we have:

$$ \|X_P(\cdot; t, x) - X_Q(\cdot; t, x)\|_{\infty, [0, \delta]} \le 2\delta L_{\Lambda,M} \|P - Q\|_{\Omega_{\delta,M}}. $$

*Proof.* By the definition of characteristic curves Eq. (2.2), for all $s \in [0, \delta]$, we have:

$$ X_P(s; t, x) - X_Q(s; t, x) = \int_t^s (\Lambda(P(\theta), X_P(\theta; t, x)) - \Lambda(Q(\theta), X_Q(\theta; t, x))) d\theta, $$
and thus by the Lipschitz property (H2) of $\Lambda$ on $[0, M]^N \times [0, 1]$:

$$ \begin{aligned} |X_P(s; t, x) - X_Q(s; t, x)| &\le \int_t^s |\Lambda(P(\theta), X_P(\theta; t, x)) - \Lambda(Q(\theta), X_Q(\theta; t, x))|d\theta \\ &\le \delta L_{\Lambda,M} \|P - Q\|_{\Omega_{\delta,M}} + \delta L_{\Lambda,M} \|X_P(\cdot; t, x) - X_Q(\cdot; t, x)\|_{\infty}. \end{aligned} $$

Taking the supremum over $s \in [0, \delta]$, we obtain:

$$ \begin{aligned} \|X_P(\cdot; t, x) - X_Q(\cdot; t, x)\|_{\infty} &\le \delta L_{\Lambda,M} \|P - Q\|_{\Omega_{\delta,M}} \\ &\quad + \delta L_{\Lambda,M} \|X_P(\cdot; t, x) - X_Q(\cdot; t, x)\|_{\infty}. \end{aligned} $$

Hence, with $\delta \le (2L_{\Lambda,M})^{-1}$, we have:

$$ \|X_P(\cdot; t, x) - X_Q(\cdot; t, x)\|_{\infty} \le 2\delta L_{\Lambda,M} \|P - Q\|_{\Omega_{\delta,M}}. \quad \square $$

LEMMA 2.4. For some $\delta > 0$ and $M > 0$ (explicitly given by Eq. (A.1) and (A.2)), the map $F$ defined by:

$$ \begin{array}{l} F: \Omega_{\delta,M} \to \Omega_{\delta,M} \\ (2.7) \qquad \qquad P \mapsto (F_1(P), \dots, F_N(P)) \end{array} $$

with for all $i \in \{1, \dots, N\}$, $F_i$ defined by Eq. (2.6), is a contraction mapping.

Proof. By regularity assumptions (H2) and (H3) on $\Lambda$ and $\mu$, and the dominated convergence Theorem, we can easily check that for all $P \in \Omega_{\delta,M}$, $F(P)$ is included in $C^0([0, \delta])^N$ (see details in Appendix A). Let us show that the image of $F$ is included in $\Omega_{\delta,M}$. For $P \in \Omega_{\delta,M}$, $i \in \{1, \dots, N\}$, $\delta > 0$ (to be chosen later) and $t \in [0, \delta]$, we have:

$$ \begin{aligned} F_i(P)(t) &\le \|\omega\|_{\infty} \int_0^t \left( r_0(u) + \sum_{i=1}^N r_i(u) P_i(u) \right) du + \|\rho_0\|_{\infty} \|\omega\|_{\infty} \\ &\le \|\omega\|_{\infty} \delta \left( \|r_0\|_{\infty, [0,T]} + M \sum_i \|r_i\|_{\infty, [0,T]} \right) + \|\rho_0\|_{\infty} \|\omega\|_{\infty}. \end{aligned} $$

Thus, for $\delta \le \frac{1}{\|\omega\|_{\infty} (\|r_0\|_{\infty, [0,T]} + M \sum_i \|r_i\|_{\infty, [0,T]})}$ and $M \ge 1 + \|\rho_0\|_{\infty} \|\omega\|_{\infty}$, $\|F_i(P)\|_{C^0([0, \delta])} \le M$ and $F(P) \in \Omega_{\delta,M}$. Now, let us show that $F$ is a contraction
mapping. Let $P, Q \in \Omega_{\delta,M}$, $i \in \{1, \dots, N\}$, and $t \in [0, \delta]$:

$$
\begin{equation}
\begin{aligned}
& \left| F_i(P)(t) - F_i(Q)(t) \right| \le \\
& \quad \int_0^t |r(u, P(u)) - r(u, Q(u))| \mathcal{F}_P(t; u, 0) \omega_i(X_P(t; u, 0)) du \\
& \quad + \int_0^t r(u, Q(u)) |\mathcal{F}_P(t; u, 0) - \mathcal{F}_Q(t; u, 0)| \omega_i(X_P(t; u, 0)) du \\
& \quad + \int_0^t r(u, Q(u)) \mathcal{F}_Q(t; u, 0) |\omega_i(X_P(t; u, 0)) - \omega_i(X_Q(t; u, 0))| du \\
& \quad + \left| \int_{1-\int_0^t \Lambda((Q(\theta), X_Q(\theta; t, 1))d\theta}^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta} \rho_0(x) \mathcal{F}_P(t; 0, x) \omega_i(X_P(t; 0, x)) dx \right| \\
& \quad + \int_0^t \rho_0(x) |\mathcal{F}_P(t; 0, x) - \mathcal{F}_Q(t; 0, x)| \omega_i(X_P(t; 0, x)) dx \\
& \quad + \int_0^t \rho_0(x) \mathcal{F}_Q(t; 0, x) |\omega_i(X_P(t; 0, x)) - \omega_i(X_Q(t; 0, x))| dx.
\end{aligned}
\tag{2.8}
\end{equation}
$$

By assumptions (H2),(H3),(H4),(H5) and (H6), all the terms on the right-hand side of Eq. (2.8) can be bounded above by a term $C(\delta) \|P - Q\|$ with $C$ a continuous function such that $C(\delta) \xrightarrow{\delta \to 0} 0$. For instance, for the fourth term on the right-hand side, using Lemma 2.3:

$$
\begin{align*}
& \left| \int_{1-\int_0^t \Lambda((Q(\theta), X_Q(\theta; t, 1))d\theta}^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta} \rho_0(x) \mathcal{F}_P(t; 0, x) \omega_i(X_P(t; 0, x)) dx \right| \\
& \le \|\rho_0\|_{\infty} \|\omega\|_{\infty} \int_0^t |\Lambda((Q(\theta), X_Q(\theta; t, 1)) - \Lambda((P(\theta), X_P(\theta; t, 1)))| d\theta \\
& \le \|\rho_0\|_{\infty} \|\omega\|_{\infty} L_{\Lambda,M} \int_0^t \|Q(\theta) - P(\theta)\|_{\mathbb{R}^N} + |X_Q(\theta; t, 1) - X_P(\theta; t, 1)| d\theta \\
& \le t \|\rho_0\|_{\infty} \|\omega\|_{\infty} L_{\Lambda,M} \left( \|P - Q\|_{\Omega_{\delta,M}} + \|X_P(\cdot; t, 1) - X_Q(\cdot; t, 1)\|_{\infty, [0, \delta]} \right) \\
& \le \delta \|\rho_0\|_{\infty} \|\omega\|_{\infty} L_{\Lambda,M} \left( \|P - Q\|_{\Omega_{\delta,M}} + 2\delta L_{\Lambda,M} \|P - Q\|_{\Omega_{\delta,M}} \right). \quad \Box
\end{align*}
$$

See Appendix A for the detailed calculations of all the terms in Eq. (2.8), which lead to the appropriate definitions of $\delta$ and $M$. Finally, by choosing a sufficiently small $\delta$ such that $C(\delta) \le \frac{1}{2} < 1$, the map $F$ is a contraction mapping on $\Omega_{\delta,M}$.

*Proof of Theorem 2.2.* From Lemma 2.4 and by the Banach fixed-point Theorem, the mapping $F$ admits a unique fixed point in $\Omega_{\delta,M}$. From that fixed point, we can construct a unique local solution to Eq. (1.1) up to time $\delta$. By iterating the process, we can extend the unique local solution on $[0, \delta_0 := \delta]$ to a unique solution on $[0, \delta_0 + \dots + \delta_n]$, where $n \in \mathbb{N}$ denotes the iteration number. In practice, at each iteration $n$, we construct a local solution $\rho^n$ of Eq. (1.1) for $t \in [0, \delta_n]$, where the initial condition is given by $\rho_0^n = \rho(\delta_0 + \dots + \delta_{n-1}, \cdot)$ and the boundary condition by $\Lambda(P^n(t))\rho^n(t, 0) = r(\delta_0 + \dots + \delta_{n-1} + t, P^n(t))$. Then we extend the solution $\rho$ by setting $\rho(\delta_{n-1} + t, \cdot) = \rho^n(t, \cdot)$ for all $t \in [0, \delta_n]$. Time $\delta_n$ is then given by the same
formula (A.2) as $\delta$, yet with $\rho(\delta_0 + \dots + \delta_{n-1}, \cdot)$ as initial data instead of $\rho_0$. To make sure that we can define a solution on $[0, T]$ for all positive $T$, we have to verify that $T$ can be reached in a finite number of iteration steps. This is the case if there is a common positive lower bound to all the $\delta'_n$s that does not depend on $n$. From condition (A.2), we can see that such a lower bound can be deduced if one obtains an upper bound on $\|\rho\|_{\infty, [0,T]}$. Let $\tau > 0$ such that the local solution $\rho$ of Eq. (1.1) is well defined up to $\tau$, thanks to the iterative construction above. For all $i \in \{1, \dots, N\}$, $P_i$ is a fixed point of Eq. (2.6). Summing all the equations, we have:

$$ 
\begin{aligned}
\sum_{i=1}^{N} P_i(t) &= \sum_{i=1}^{N} \int_{0}^{t} r(u, P(u)) \mathcal{F}_P(t; u, 0) \omega_i(X_P(t; u, 0)) du \\
&\quad + \sum_{i=1}^{N} \int_{0}^{1-\int_{0}^{t} \Lambda((P(\theta), X_P(\theta; t, 1))) d\theta} \rho_0(x) \mathcal{F}_P(t; 0, x) \omega_i(X_P(t; 0, x)) dx \\
&\leq B + C \int_{0}^{t} \sum_{i=1}^{N} P_i(u) du,
\end{aligned}
$$

with $B = N \|\omega\|_{\infty} (\|\rho_0\|_{\infty} + T \|r_0\|_{\infty})$ and $C = N \|\omega\|_{\infty} \sup_i \|r_i\|_{\infty}$. Then by Grönwall's Lemma, for all $t \in [0, \tau]$ and $i \in \{1, \dots, N\}$,

$$ (2.9) \qquad P_i(t) \leq \sum_{i=1}^{N} P_i(t) \leq Be^{CT} \leq M, $$

where $M$ is given by (A.1). By the quasi-explicit formula Eq. (2.5), and using (2.3), for all $t \in [0, T]$:

$$ 
\begin{cases} 
|\rho(t, x)| \leq \|\rho_0\|_{\infty} e^{-\int_0^t \Lambda_x(P(\theta), X_P(\theta; t, x))} , & X_{0,P(t)} \leq x \leq 1, \\ 
|\rho(t, x)| \leq \frac{r(\alpha_P(t, x), P(\alpha_P(t, x)))}{\Lambda(P(\alpha_P(t, x)), 0)} e^{-\int_{\alpha_P(t, x)}^t \Lambda_x(P(\theta), X_P(\theta; t, x)) d\theta} , & 0 \leq x \leq X_{0,P(t)}. 
\end{cases}
$$

Finally, using the upper bound (2.9), we can express a bound on $\rho$ that only depends on $T$:

$$ \|\rho\|_{\infty, [0,1] \times [0,\tau]} \le \exp \left( T \|\Lambda_x\|_{\infty, [0,M]^N \times [0,1]} \right) \max \left( \|\rho_0\|_{\infty}, \frac{\|r\|_{\infty, [0,T] \times [0,M]^N}}{\min_{[0,M]^N \times [0,1]} \Lambda} \right). $$

In particular, for all $n \in \mathbb{N}$, $\rho_0^n$ is bounded uniformly in $n$ as long as $\sum_{i \le n} \delta_i < T$. $\square$

**2.3. Stationary solutions.** From now on, unless otherwise stated, we will make the following assumption:

(H7) The recruitment rate $r$ is autonomous, i.e. $r(t, P) = r(P)$.

Stationary solutions to Eq. (1.1) are the functions $\bar{\rho} \in C^1([0, 1], \mathbb{R}^+)$ verifying the following problem:

$$ (2.10) \qquad \left\{ \begin{array}{l} \bar{\rho}'(x) = -\frac{\mu(\bar{P}, x) + \Lambda_x(\bar{P}, x)}{\Lambda(\bar{P}, x)} \bar{\rho}(x), \quad x \in [0, 1], \\ \bar{\rho}(0) = \frac{r(\bar{P})}{\Lambda(\bar{P}, 0)}, \end{array} \right. $$
where

$$ \bar{P} = (\bar{P}_i)_{1\le i \le N}, \quad \bar{P}_i = \int_0^1 \omega_i(x) \bar{p}(x) dx. $$

In the same spirit as for the dynamic problem, we can limit the study of the differential Eq. (2.10) to an equivalent integral equation on $\bar{P}$. This is the statement of the following proposition, which directly comes from the quasi-explicit formula verified by the solution to Eq. (2.10) and the regularity assumptions (H2), (H3), (H4).

**PROPOSITION 2.5.** *Stationary solutions $\bar{p}$ of Eq. (1.1) are defined by:*

$$ \bar{p}(x) = \frac{r(\bar{P})}{\Lambda(\bar{P}, x)} e^{-\int_0^x \frac{\mu(\bar{P}, y)}{\Lambda(\bar{P}, y)} dy} \quad \text{on } (0, 1), $$

where $\bar{P}$ is a fixed point of $F: (\mathbb{R}^+)^N \mapsto (\mathbb{R}^+)^N$ defined by:

$$ (2.11) \qquad F_i(P) = \int_0^1 \omega_i(x) \frac{r(P)}{\Lambda(P, x)} e^{-\int_0^x \frac{\mu(P, y)}{\Lambda(P, y)} dy} dx, \quad i \in \{1, \dots, N\}. $$

*Stationary solutions to Eq. (1.1) are in one-to-one correspondence with fixed points of F.*

The number of stationary solutions can vary from zero to infinity, depending on the formulation of the growth, death, and recruitment rates. Nevertheless, if we assume that the recruitment-to-growth ratio is bounded, then we can show the existence of at least one stationary solution.

**PROPOSITION 2.6.** *Assuming that:*

$$ (2.12) \qquad \exists M, \forall P \in (\mathbb{R}^+)^N, \forall x \in [0, 1], \frac{r(P)}{\Lambda(P, x)} \le M, $$

then Eq. (2.10) has at least one solution.

*Proof.* Due to Proposition 2.5, we only have to show that, under the assumption given by Eq. (2.12), $F$ defined by Eq. (2.11) has at least one fixed point. First, note that, thanks to the regularity assumptions (H2), (H3), and (H4) on $\Lambda$, $\mu$ and $r$, $F$ is continuous. Let $C := M \times \max_{1 \le i \le N} \|\omega_i\|_{\infty}$, then for all $i \in \{1, \dots, N\}$,

$$ F_i(P) = \int_0^1 \omega_i(x) \frac{r(P)}{\Lambda(P, x)} e^{-\int_0^x \frac{\mu(P, y)}{\Lambda(P, y)} dy} dx \le \|\omega_i\|_{\infty} \int_0^1 \frac{r(P)}{\Lambda(P, x)} dx \le C. $$

The image of $F$ is in $[0, C]^N$, thus the restriction of $F$ to $[0, C]^N$ is a continuous application from a compact convex set to itself. By Brouwer's fixed-point Theorem, $F$ has at least one fixed point in $[0, C]^N$. □

To discuss stationary solutions in more details, the model needs to be refined. In the following section, under the key assumption of a separable growth rate and in the absence of death, this refinement will enable us to better characterize stationary solutions and study their stability.
**3. Asymptotic analysis with a separable growth rate.** In this section, we make three additional assumptions:

(H8) The growth rate is separable: $\Lambda(P, x) = f(x)g(P)$, with $f$ and $g$ strictly positive functions.

(H9) The death rate is negligible: $\mu(P, x) \equiv 0$.

(H10) The ratio $b := \frac{r}{g}$ is bounded on $\mathbb{R}_+^N$, and $g$ admits a strictly positive infimum on $\mathbb{R}_+^N$. Furthermore, $g$ and $r$ are continuously differentiable and globally Lipschitz continuous on $\mathbb{R}_+^N$.

Then, under (H7)-(H10), Eq. (1.1) becomes:

$$ (3.1) \qquad \begin{cases} \partial_t \rho + g(P(t)) \partial_x (f(x)\rho(t,x)) = 0, & x \in [0, 1], t \ge 0 \\ f(0)g(P(t))\rho(t,0) = r(P(t)), & t \ge 0 \\ \rho(0,x) = \rho_0(x), & x \in [0, 1]. \end{cases} $$

The key additional assumption is (H8), as it allows us to perform a time-scaling transformation. Consequently, studying Eq. (3.1) will become equivalent to studying an equation with a growth rate independent of $P$ (actually even a constant), namely a semilinear one.

Assumption (H9) aligns well with our main biological application. Indeed, in many fish species, oocyte death is not readily apparent under normal physiological conditions. While not strictly needed, from a mathematical perspective, assumption (H9) simplifies calculations, particularly the characterization of stationary solutions (see Proposition 2.5 against Proposition 3.1) and their stability. Indeed without (H9), Proposition 3.6 is still true but with a different characteristic equation which is a little bit harder to analyze (see for example [10] for a similar model with death and growth rates that do not depend on $x$).

Assumption (H10) ensures the existence of at least one stationary solution (by Proposition 2.6), and guarantees that function $\frac{r}{g}$ is continuously differentiable and globally Lipschitz continuous on $\mathbb{R}_+^N$. Moreover, this assumption will be helpful for studying the stability of stationary solutions.

Let us note that with the additional assumptions (H8)-(H10), it can be shown that a unique weak solution to Eq. (3.1) in $C^0([0, T]; L^1(0, 1))$ exists for any $L^1$ initial condition (instead of $L^\infty$ in Theorem 2.2). We derive a sketch of the proof adapted to this case in Appendix B.

As a direct application of Proposition 2.5 under assumptions (H8) and (H9), we can provide a characterization of stationary solutions to Eq. (3.1).

**PROPOSITION 3.1.** *Under assumptions (H1)-(H10), the stationary solutions $\bar{\rho}$ of Eq. (3.1) are defined by:*

$$ \bar{\rho}(x) = \frac{b(\bar{P})}{f(x)}, \quad x \in [0, 1], \quad \text{where } \bar{P} \text{ verifies } \bar{P} = b(\bar{P})T \text{ in } \mathbb{R}_+^N, $$

$$ \textit{and with } T := \int_0^1 \frac{\omega}{f} \in \mathbb{R}_+^N. $$

We emphasize that $\bar{P} = b(\bar{P})T$ can be rewritten $\bar{P} = r(\bar{P}) \int_0^1 \frac{\omega(x)}{f(x)g(\bar{P})} dx$, where $r(\bar{P})$ is the number of new individuals per unit of time at equilibrium, and $\int_0^1 \frac{\omega(x)}{f(x)g(\bar{P})} dx$ is the weighted lifetime of an individual. At equilibrium, it is natural that the weighted population average $\bar{P}$ is equal to the number of new individuals per unit of time
multiplied by their weighted lifetime. Let us also point out that all the stationary solutions to Eq. (3.1) have the same shape as $\frac{1}{f}$.

**3.1. Time and space scaling transformation.** We denote by

$\Gamma(x) := \int_{0}^{x} \frac{1}{f(y)} dy$ the function that measures the time required for an individual to grow from initial size 0 to size $x$ in the absence of interaction terms. We further denote $\tau := \Gamma(1)$ as the entire lifetime of an individual in the absence of interaction terms. Then the change of variables:

$$ (3.2) \qquad s = \int_{0}^{t} g(P(\theta))d\theta $$

$$ (3.3) \qquad y = \Gamma(x) $$

lead to the following problem:

$$ (3.4) \qquad \left\{ \begin{array}{l} \partial_s v(s, y) + \partial_y v(s, y) + m(y)v(s, y) = 0, \quad y \in (0, \tau), \quad s \ge 0 \\ v(s=0, y) = v_0(y), \quad y \in (0, \tau) \\ v(s, 0) = B(v(s, \cdot)), \quad s \ge 0, \end{array} \right. $$

where $m := f' \circ \Gamma^{-1}$, $\alpha := (\omega \circ \Gamma^{-1})(f \circ \Gamma^{-1})$ and for all $\phi \in L^1(0, \tau)$, $B(\phi) := \frac{1}{f(0)}b(\int_{0}^{\tau} \alpha\phi)$. We use the same notion of weak solution to Eq. (3.4) as that provided in Definition 2.1. By assumptions (H2), (H6), and (H10) we note that $\alpha$, $m$ and $B$ satisfy:

(H11) For all $i \in [1; N]$, $\alpha_i : [0, \tau] \mapsto \mathbb{R}$ is nonnegative and Lipschitz continuous.

(H12) The function $m : [0, \tau] \mapsto \mathbb{R}$ is continuous.

(H13) The map $B : L^1(0, \tau) \mapsto \mathbb{R}_+^N$ is nonnegative, continuously differentiable and Lipschitz continuous.

Let $T \ge 0$. For all $\rho \in C([0, T], L^1(0, 1))$, we define

$$ (3.5) \qquad \begin{array}{l} s_{\rho} : [0, T] \to \mathbb{R} \\ t \mapsto \int_{0}^{t} g(P(\theta))d\theta \end{array} $$

with $P(\theta) = \int_{0}^{1} \omega(x)\rho(\theta, x)dx$.

Let $S \ge 0$. For all $v \in C([0, S], L^1(0, \tau))$, we define

$$ (3.6) \qquad \begin{array}{l} t_v : [0, S] \to \mathbb{R} \\ s \mapsto \int_{0}^{s} \frac{1}{g(V(\theta))} d\theta \end{array} $$

with $V(\theta) = \int_{0}^{\tau} \alpha(x)v(\theta, x)dx$.

The following three propositions establish the link between Eq. (3.1) and Eq. (3.4).

**PROPOSITION 3.2.** Under assumptions (H1)-(H10), if $\rho$ is the weak solution to Eq. (3.1) on $[0, T] \times [0, 1]$, then the only $v \in C^0([0, s_\rho(T)], L^1(0, \tau))$ that verifies $v(s_\rho(t), \Gamma(x)) = \rho(t, x)$ for all $(t, x) \in [0, T] \times [0, 1]$ is the weak solution to Eq. (3.4) on $[0, s_\rho(T)] \times [0, \tau]$ with initial condition $v_0 = \rho_0 \circ \Gamma^{-1}$.

Conversely, if $v$ is the weak solution to Eq. (3.4) on $[0, S] \times [0, \tau]$, then the only $\rho \in C^0([0, t_v(S)], L^1(0, 1))$ that verifies $\rho(t_v(s), \Gamma^{-1}(y)) = v(s, y)$ for all $(s, y) \in [0, S] \times [0, \tau]$ is the weak solution to Eq. (3.1) on $[0, t_v(S)] \times [0, 1]$ with initial condition $\rho_0 = v_0 \circ \Gamma$.
The proof of Proposition 3.2 is given in Appendix C and relies on the change of variables Eq. (3.2)-(3.3) and assumption (H2) on the growth rate.

Since we are interested in the asymptotic behavior of the solutions, we want to guarantee that stationary solutions to Eq. (3.1) and Eq. (3.4) are mapped, and that stability is preserved. We recall the standard notion of local and asymptotic stability in Definition C.1.

**PROPOSITION 3.3.** Under assumptions (H1)-(H10), $\bar{\rho}$ is a stationary solution to Eq. (3.1) if and only if $\bar{v} := \bar{\rho} \circ \Gamma^{-1}$ is a stationary solution to Eq. (3.4). In such a case, $\bar{P} = \int_0^1 \omega(x)\bar{\rho}(x)dx = \int_0^\tau \alpha(x)\bar{v}(x)dx = \bar{V}$. Furthermore, $\bar{\rho}$ is locally $L^1$-stable (resp. locally $L^1$-asymptotically stable) if and only if $\bar{v}$ is locally $L^1$-stable (resp. locally $L^1$-asymptotically stable).

*Proof.* The first assertion is immediate by Proposition 3.2. Let us show that the second assertion holds, i.e. that stability is preserved. Let us assume that $\bar{\rho}$ is a locally $L^1$-stable stationary solution to Eq. (3.1). Let $\epsilon > 0$. By definition, there exists $\delta_\rho > 0$ such that
(3.7)
$$ \forall \rho_0 \in L^1(0,1), \quad (\|\bar{\rho} - \rho_0\|_{L^1(0,1)} \le \delta_\rho) \Rightarrow \left(\forall t \ge 0, \|\bar{\rho} - \rho(t, \cdot)\|_{L^1(0,1)} \le \epsilon \min_{[0,1]} f\right). $$

Let $v_0 \in L^1(0, \tau)$ such that $\|\bar{v} - v_0\|_{L^1(0,\tau)} \le \frac{\delta_\rho}{\max_{[0,\tau]} f \circ \Gamma^{-1}}$. Define $\rho_0 = v_0 \circ \Gamma$, then

$$ \|\bar{\rho} - \rho_0\|_{L^1(0,1)} = \|\bar{v} \circ \Gamma - v_0 \circ \Gamma\|_{L^1(0,1)} \le \max_{[0,\tau]} f \circ \Gamma^{-1} \|\bar{v} - v_0\|_{L^1(0,\tau)} \le \delta_\rho. $$

Therefore, for $s \ge 0$, using Eq. (3.7), we deduce:

$$ \begin{gathered} \|\bar{v} - v(s)\|_{L^1(0,\tau)} = \|\bar{\rho} \circ \Gamma^{-1} - \rho(t_v(s), \Gamma^{-1}(\cdot))\|_{L^1(0,\tau)} \\ (3.8) \qquad \le \frac{1}{\min_{[0,1]} f} \|\bar{\rho} - \rho(t_v(s), \cdot)\|_{L^1(0,1)} \le \epsilon, \end{gathered} $$

and $\bar{v}$ is locally $L^1$-stable. In such a case, the map $\theta \mapsto V(\theta)$ is bounded on $[0, +\infty)$, and by continuity, $g \circ V$ has a strictly positive maximum, so that $t_v(s) \ge Cs$ for some $C > 0$ (see Eq. (3.6)). Thus, $t_v(s) \xrightarrow{s \to +\infty} +\infty$, and the same calculus as in Eq. (3.8) shows that:

$$ \|\bar{v} - v(s)\|_{L^1(0,\tau)} \xrightarrow{s \to +\infty} 0, $$

whenever $\|\bar{\rho} - \rho(t, \cdot)\|_{L^1} \xrightarrow{t \to +\infty} 0$.

The reciprocal works in the same way.

□

By directly applying Proposition 3.2, we can also state the following.

**PROPOSITION 3.4.** If $\rho$ is a T-periodic solution to Eq. (3.1), then the associated solution $v$ of Eq. (3.4) given by Proposition 3.2 is $s_\rho(T)$-periodic. Conversely, if $v$ is a S-periodic solution to Eq. (3.4), then the associated solution $\rho$ of Eq. (3.1) given by Proposition 3.2 is $t_v(S)$-periodic.

Furthermore, let us point out that with the same arguments as those used in the proof of Proposition 3.3 for the stability part, we expect that orbital $L^1$-stability of periodic solutions is preserved by the transformation.
Propositions 3.3 and 3.4 allow us to study the stability of stationary solutions and the existence of periodic solutions to Eq. (3.1) by analyzing these properties in Eq. (3.4), which is simpler than the original semilinear equation. This will be the matter of the next section.

**3.2. Asymptotic analysis via local linearization.** Eq. (3.4) can be rewritten as a non-densely defined semilinear abstract Cauchy problem, in the framework of [37]. We consider the Banach space $X := \mathbb{R} \times L^1(0, \tau)$ endowed with the norm $\|(a, \phi)\|_X := |a| + \|\phi\|_{L^1(0, \tau)}$.
Let us define the following linear operator $A: D(A) \subset X \to X$ by

$$ (3.9) \qquad A\begin{pmatrix} 0 \\ \phi \end{pmatrix} = \begin{pmatrix} -\phi(0) \\ -\partial_y\phi - m\phi \end{pmatrix}, \quad D(A) = \{0\} \times W^{1,1}(0, \tau) $$

and the map $G: \overline{D(A)} \to X$ by

$$ G\begin{pmatrix} 0 \\ \phi \end{pmatrix} = \begin{pmatrix} B(\phi) \\ 0 \end{pmatrix}. $$

To simplify the notation, set $v(t) = \begin{pmatrix} 0 \\ v(t, \cdot) \end{pmatrix}$. Equation (3.4) can then be rewritten as the following abstract Cauchy problem:

$$ (3.10) \qquad \frac{dv(t)}{dt} = Av(t) + G(v(t)), \quad t \ge 0, \quad v(0) = \begin{pmatrix} 0 \\ v_0 \end{pmatrix} \in \overline{D(A)}. $$

Thanks to (H13), $G$ is continuously differentiable and Lipschitz continuous.

The linearized problem around stationary solution $\begin{pmatrix} 0 \\ \bar{v} \end{pmatrix}$, abbreviated $\bar{v}$ for brevity, is then:

$$ \frac{dw}{dt} = Aw(t) + DG(\bar{v})w(t), \quad t \ge 0, \quad w(0) = w_0 \in \overline{D(A)}, $$

where $DG(\bar{v}) \in \mathcal{L}(\overline{D(A)}, X)$ is the Fréchet derivative of $G$ at $\bar{v}$, and is given explicitly as:

$$ (3.11) \qquad DG(\bar{v})\begin{pmatrix} 0 \\ \phi \end{pmatrix} = \begin{pmatrix} \frac{1}{f(0)} \nabla b(\bar{V}) \cdot \int_0^\tau \alpha \phi \\ 0 \end{pmatrix}, \quad \forall \phi \in L^1(0, \tau), $$

with $\bar{V} = \int_0^\tau \alpha \bar{v}$. Note that $A$ is not-densely defined, since $\overline{D(A)} = \{0\} \times L^1(0, \tau) \neq X$. To apply the results of [37] on the stability of equilibria and Hopf bifurcation Theorem, we need to gather stability estimates on $A$, and further need some control on the perturbation $DG(\bar{v})$.

First of all, let us recall some usual notations from spectral theory (see Chap. 4 in [18] for further details). Let $L: D(L) \subset X \to X$ a closed linear operator. We denote by $\rho(L) := \{\lambda \in \mathbb{C} | \lambda - L$ is bijective\}$ the resolvent set of $L$, $\sigma(L) := \mathbb{C} \setminus \rho(L)$ its spectrum, and $s(L) := \{\sup \Re(\lambda) | \lambda \in \sigma(L)\}$ its spectral bound. The operator norm of $L$ is defined by $\|L\|_{L(X)} := \sup_{x \in D(L) \setminus \{0\}} \frac{\|Lx\|_X}{\|x\|_X}$. Let $\Omega \subset X$, we note $L_\Omega$ the part of $L$ in $\Omega$, defined by $L_\Omega x := Lx, \forall x \in D(L_\Omega) := \{x \in D(L) : Lx \in \Omega\}$.
PROPOSITION 3.5. Let assumptions (H11)-(H13) hold. A, defined by Eq. (3.9), is a closed Hille-Yosida operator,

$$(-\underline{m}, +\infty) \subset \rho(A)$$

and

$$(3.12) \qquad \|(\lambda \mathbb{I} - A)^{-n}\|_{\mathcal{L}(X)} \le \frac{M}{(\lambda + \underline{m})^n}, \quad \forall \lambda > -\underline{m}, \forall n \ge 1$$

where $\underline{m} := \inf_{y \in [0, \tau]} m(y) > -\infty$.

Let $\bar{v}$ a stationary solution to Eq. (3.10), and $DG(\bar{v})$ defined by Eq. (3.11). $DG(\bar{v})$ is a compact operator from $\overline{D(A)}$ to $X$. Furthermore, the essential growth bound of $(A + DG(\bar{v}))_{\overline{D(A)}}$ is strictly negative, the spectral bound of $(A + DG(\bar{v}))_{\overline{D(A)}}$ equals the growth bound of its associated semigroup, and the spectrum of $A + DG(\bar{v})$ is only made up of eigenvalues, which can be determined by the zeros of a characteristic function, i.e.:

$$\sigma(A + DG(\bar{v})) = \sigma_p(A + DG(\bar{v})) = \{\lambda \in \mathbb{C} : \Delta(\lambda) = 0\},$$

with

$$(3.13) \qquad \Delta(\lambda) := 1 - \nabla b(\bar{V}) \cdot \int_0^1 \frac{\omega(x)}{f(x)} e^{-\lambda \Gamma(x)} dx,$$

and $\cdot$ denotes the usual dot product of $\mathbb{R}^N$. Moreover, the characteristic function $\Delta$ defined by Eq. (3.13) is holomorphic on $\mathbb{C}$, and for all $n \in \mathbb{N}^*$,

$$\frac{d^n \Delta(\lambda)}{d\lambda^n} = (-1)^{n+1} \nabla b(\bar{V}) \cdot \int_0^1 \frac{\omega(x)}{f(x)} \Gamma(x)^n e^{-\lambda \Gamma(x)} dx.$$

Finally, the roots of $\Delta$ are isolated, and $\lambda$ is a simple isolated eigenvalue of $A+DG(\bar{v})$ if and only if $\frac{d\Delta(\lambda)}{d\lambda} \ne 0$.

Proof. The properties on A are rather standard and the proof leading to Eq. (3.12) follows from a direct calculation of the resolvent, which involves a first-order ordinary differential equation, and is given in Appendix D. As $\alpha$ is bounded, the image of the unit ball of $\overline{D(A)}$ by $DG$ is included in $[-M, M] \times \{0, L_1(0, \tau)\}$ for some $M$, which is compact in $X$. As A is a Hille-Yosida operator, and $DG(\bar{v})$ is a compact operator, then for all $t \ge 0$, $(DG(\bar{v})) \circ T_{A_0}(t)$ is compact (as the composition of a compact operator with a bounded linear operator is compact). We can thus apply Theorem 1.2. in [17] to get that the essential growth bound of $(A + DG(\bar{v}))_{\overline{D(A)}}$ is less than or equal to the essential growth bound of $A_0 := A_{\overline{D(A)}}$, the part of A in $\overline{D(A)}$. Using a corollary of the Hille-Yosida Theorem (see Lemma 3.4.2. in [37]), $A_0$ generates a $C_0$-semigroup $T_{A_0}$ of bounded linear operators on $\overline{D(A)}$, with $\|T_{A_0}(t)\|_{\mathcal{L}(\overline{D(A)})} \le e^{-mt}$, $\forall t \ge 0$. We can even specify $T_{A_0}$ directly from the definition of $A_0$:

$$T_{A_0}(t) \begin{pmatrix} 0 \\ \phi \end{pmatrix} = \begin{pmatrix} 0 \\ \hat{T}_{A_0}(t)\phi \end{pmatrix}$$

where

$$(3.14) \qquad \hat{T}_{A_0}(t)\phi(y) = \begin{cases} e^{-\int_{y-t}^y m(x)dx} \phi(y-t), & \text{if } 0 \le t \le y \le \tau \\ 0 & \text{otherwise.} \end{cases}$$
We deduce from Eq. (3.14) that for every $t > \tau$, $T_{A_0}(t) = 0$. Hence $T_{A_0}$ is a nilpotent $C_0$-semigroup and the growth bound of $T_{A_0}$ is $\omega_0(A_0) = -\infty$. By [18] Proposition IV.2.10., the essential growth bound of $A_0$ is less than the growth bound of $T_{A_0}$, so that the essential growth bound is equal to $-\infty < 0$.

Furthermore, it is well known (corollary IV.2.11 in [18]) that the growth bound of a $C_0$-semigroup is the maximum between the essential growth bound and the spectral bound of its generator. Finally, we determine the spectrum of $A + DG(\bar{v})$ by an explicit calculation of its resolvent, given in Appendix D. Properties of $\Delta$ follows from similar arguments as those used to derive Lemma 8.3.6. in [37]. □

We can now state the main result of this section, which enables us to study the stability of stationary solutions to Eq. (3.1) by analyzing the roots of a characteristic function.

**PROPOSITION 3.6.** Let $\bar{\rho}$ a stationary solution to Eq. (3.1). Under assumptions (H10)-(H13), if the complex roots of the characteristic function Eq. (3.13) are all strictly negative, then $\bar{\rho}$ is locally asymptotically stable. If there exists a complex root with a positive real part, then $\bar{\rho}$ is unstable.

*Proof.* The proof is a direct application to Eq. (3.10) of Propositions 5.7.1. and 5.7.4. in [37] together with Proposition 3.5. Then Proposition 3.3 transfers the result to the stationary solutions to Eq. (3.1). □

Under certain conditions, the location of the roots of the characteristic function Eq. (3.13) is easier to determine. These results will be summarized in the following corollaries.

**COROLLARY 3.7.** Let $\bar{\rho}$ a stationary solution to Eq. (3.1). Under assumptions (H10)-(H13), and if $b$ is locally increasing near $\bar{P}$, i.e.,

$$ (3.15) \qquad \nabla b(\bar{P}) \ge 0, $$

then:

* $\bar{\rho}$ is locally asymptotically stable if $\mathcal{T} \cdot \nabla b(\bar{P}) < 1$,

* $\bar{\rho}$ is unstable if $\mathcal{T} \cdot \nabla b(\bar{P}) > 1$.

*Proof.* From Proposition 3.3, let's consider $\bar{v}$ the stationary solution associated with Eq. (3.4). We show in Appendix D that if $\nabla b(\bar{v}) \ge 0$ (recall that $\bar{P} = \bar{V}$), the $C_0$-semigroup generated by $A + DG(\bar{v})$ is positive. The positivity of the semigroup implies that the spectral bound $s(A + DG(\bar{v}))$ is real by direct application of Theorem VI.1.10 in [18]. Hence the study of the characteristic function defined by Eq. (3.13) is greatly simplified and proceeds as follows. The restriction to $\mathbb{R}$ of $\Delta$, $\lambda \in \mathbb{R} \mapsto \Delta(\lambda) \in \mathbb{R}$ is a strictly increasing function with $\lim_{\lambda \to \infty} \Delta(\lambda) = -\infty$, $\lim_{\lambda \to \infty} \Delta(\lambda) = 1$, so that it admits a unique root, which is strictly positive if and only if $\Delta(0) = 1 - \mathcal{T} \cdot \nabla b(\bar{P}) < 0$. Then by applying Proposition 3.6, we get the result. □

**COROLLARY 3.8.** Let $\bar{\rho}$ a stationary solution to Eq. (3.1). Under assumptions (H10)-(H13), and if $\omega(x) \equiv \omega \in \mathbb{R}^N$ is constant, then:

* $\bar{\rho}$ is locally asymptotically stable if $\mathcal{T} \cdot \nabla b(\bar{P}) < 1$,

* $\bar{\rho}$ is unstable if $\mathcal{T} \cdot \nabla b(\bar{P}) > 1$, and the characteristic function admits a positive real root.

Also, if $\mathcal{T} \cdot \nabla b(\bar{P}) = 1$, then the only root of $\Delta(\lambda)$ is 0.
*Proof.* From Eq. (3.13), we can establish that $\Delta(\lambda) = 0$ if and only if

$$ \left( \lambda = 0 \text{ and } \nabla b(\bar{P}) \cdot \omega = \frac{1}{\tau} \right) \text{ or } \left( \lambda \neq 0 \text{ and } \frac{1 - e^{-\lambda \tau}}{\lambda} \nabla b(\bar{P}) \cdot \omega = 1 \right). $$

Then, we can use Theorem 2 in [26] to handle with the case $\lambda \neq 0$, which is a well-known transcendental equation in the field of delay differential equations. □

The following result is derived from Corollaries 3.7 and 3.8:

**COROLLARY 3.9.** Under assumptions (H10)-(H13), and assuming that $r(0_{\mathbb{R}^N}) > 0$ (to ensure the existence of at least one non-trivial stationary solution), and that Eq. (3.1) has a finite number of stationary solutions, then:

* If $\omega(x) \equiv \omega \in \mathbb{R}^N$ is constant or if $b$ is non decreasing with respect to all variables, then model Eq. (3.1) has $k$ unstable hyperbolic stationary solutions, and $k + 1$ locally asymptotically stable hyperbolic stationary solutions, with $k \ge 0$ an integer which depends on $b$.

* If $b$ is non increasing with respect to all variables, then there exists a unique stationary solution.

*Proof.* Let us define $q(X) = \mathcal{T}_1 b(X, \frac{T_2}{T_1} X, \dots, \frac{T_N}{T_1} X) - X$ for $X \ge 0$. Then from Proposition 3.1, $\bar{P} = (\bar{P}_1, \bar{P}_2, \dots, \bar{P}_N) = (\bar{P}_1, \frac{T_2}{T_1} \bar{P}_1, \dots, \frac{T_N}{T_1} \bar{P}_1)$ defines a unique stationary solution $\bar{\rho}(x) = \frac{b(\bar{P})}{f(x)}$ if and only if $q(\bar{P}_1) = 0$. Moreover, by composition $q'(\bar{P}_1) = \mathcal{T} \cdot \nabla b(\bar{P}) - 1$. Firstly, if $\omega(x) \equiv \omega \in \mathbb{R}^N$ is constant and/or if $b$ is non decreasing with respect to all variables, the number of stationary solutions is given by the number of roots of $q$ in $\mathbb{R}^+$, and by Corollary 3.7 and 3.8, the stability is given by the sign of the derivative at the roots. As $b$ is bounded, the function $q$ is negative when $X$ is large enough. Moreover, as $r(0_{\mathbb{R}^N}) > 0$, then $b(0_{\mathbb{R}^N}) > 0$ so $q(0) > 0$. Thus, $q$ crosses the abscissa line from the positive half plane to the negative half plane one time more than in the opposite direction. The first result follows. If $b$ is non increasing with respect to all variables, then $q$ is decreasing on $\mathbb{R}^+$, and as $q(0) > 0$ and $q(X) < 0$ when $X$ is large enough, $q$ has a unique root, so there is a unique stationary solution.

To end with this section, let us note that when $\omega(x) \equiv \omega \in \mathbb{R}^N$ is constant and/or if $b$ is non decreasing with respect to all variables, we do not expect oscillatory solutions. Indeed, when $b$ is non decreasing, the leading eigenvalue of the linearized operator is real (see Proof of Corollary 3.7), and, when $\omega$ is constant, transversality conditions leading to Hopf bifurcation cannot be satisfied (see the end of Corollary 3.8). On the other side, when $\omega$ is not constant and $b$ is non increasing with respect to all variables, we can expect oscillatory solutions. This behavior will be shown off in the next section.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Main steps of fish oogenesis with feedback loops. The color codes delineates the steps for which separate models are designed. The horizontal arrows correspond to the growth rate, the circle arrows correspond to the recruitment into each step. The + (resp. −) symbols discriminate positive against negative feedback loops mediated by different hormonal signals.</figcaption></figure>

**4. Applications and numerical illustrations.** In this section, we detail the application to fish oogenesis. We take advantage of the sequential progress of oocyte maturation to consider three separate (sub)models illustrating three possible configurations of the model given by Eq. (3.1), namely a double positive or negative feedback loop and a negative-positive feedback loop.

Fish (follicle-enclosed) oogenesis can indeed be divided into three main steps, the pre-vitellogenic phase (small-sized oocytes), the vitellogenic phase (medium-sized oocytes) and the post-vitellogenic phase (large-sized oocytes). These phases are named after the vitellogenesis process, during which the oocyte accumulates the lipoprotein vitellogenin that will later feed the embryo. Each phase is further characterized by the secretion of a prominent hormonal signal, respectively AMH (Anti-Müllerian Hormone), Estradiol, and the progestagen DHP (17α,20βDiHydroxy-4-Pregnen-3-one), which in turn either promotes or inhibits the recruitment into and growth within the corresponding phase, as illustrated in Figure 1. AMH acts locally on the ovarian level, while Estradiol and DHP are involved in complex endocrine loops involving the hypothalamus and pituitary gland, and also the liver through the Estradiol-Vitellogenin feedback loop. In the model formulation, we do not differentiate the AMH paracrine effect from the Estradiol/DHP endocrine effects, and all hormonal controls are computed directly as weighted integrals of the density.

In this section, we can apply the theoretical results from Sec. 3 and exhibit oscillatory and bistable solutions in a biologically-relevant context.

We consider the case of a single interaction term $P(t) = \int_0^1 \omega_i(x)\rho(t,x)dx$. In such a case, we can make the following comments on the model behavior. First, if the recruitment $r$ decreases with $P$ (negative feedback) and $g$ increases with $P$ (positive feedback), then $b = r/g$ decreases with $P$. From Corollary 3.9, there is a unique stationary solution, whose stability remains to be determined. This situation is investigated in the first application in subsection 4.1, which exemplifies the occurrence of a Hopf bifurcation. Conversely, when $b$ increases with $P$, there may be several stationary solutions as stated in Corollary 3.9. Following the proof of Corollary 3.9, we can provide a graphical interpretation of a multi-stationary behavior, as presented in the second application in subsection 4.2.

**4.1. A model of the spawning cycle (postvitello).** In teleost fish there is a variety of reproductive strategies. Some species (typically the salmonids) spawn only once (and die thereafter) and the oocyte population is synchronous at any stage.
Most other groups exhibit spawning cycles related to an asynchronous progress along oogenesis [46]. In the latter case, the most mature follicles have fast dynamics and the post-vitellogenic phase is completed within a spawning cycle.

For this application, the lower size boundary ($x = 0$) in Eq. (3.1) corresponds to the size of newly recruited large-sized (postvitellogenic) oocytes, and the greater size boundary $x = 1$ to the size of large-sized oocytes ready to be spawned.

From available data on DHP secretion pattern before spawning (deduced from the corresponding expression of the enzyme responsible for DHP secretion) [50], we model roughly $\omega(x)$ as a triangular function, with $\omega(0) = 0$, $\omega(0.5) = M > 0$, $\omega(1) = 0$:

$$ (4.1) \qquad \omega(x) = 2Mx, \quad x \in [0, 0.5], \quad \omega(x) = 2M(1-x), \quad x \in [0.5, 1]. $$

DHP promotes the growth of large-sized oocytes, which leads to a positive control on the growth rate $g$. Furthermore, through an endocrine feedback loop involving the central components of the reproductive axis (the hypothalamus and pituitary gland), DHP prevents the recruitment of new large-sized oocytes. We model these controls as Hill functions as follows:

$$ (4.2) \qquad \Lambda(x, P) = f(x)g(P) = f(x) \frac{1+c_{\lambda}P^n}{1+P^n}, \quad r(P_r) = \frac{c_r}{1+P^n} $$

with $c_{\lambda} > 1$, $c_r > 0$, $n > 0$. Parameter $c_r$ represents the basal rate of newly recruited large-sized oocyte, and $c_{\lambda}$ is the maximal growth rate of large-sized oocytes. From Eq. (4.2), the ratio $b := \frac{r}{g}$ is decreasing with $P$:

$$ (4.3) \qquad b(P) = \frac{c_r}{1+c_{\lambda}P^n}. $$

In what follows, we will study the asymptotic behavior (stationary solution existence and stability, Hopf bifurcation) of the model, with respect to parameter $c_r$.

**4.1.1. Bifurcation analysis.** By direct application of Proposition 3.1, we can derive a characterization of the model stationary solution.

**PROPOSITION 4.1.** For all $c_{\lambda} > 1$, $c_r > 0$ and $n > 0$, the spawning cycle model admits a unique stationary solution $\bar{\rho}(x) = \frac{\bar{P}}{\bar{T}f(x)}$ on $[0, 1]$, where $\bar{P}$ is the unique strictly positive solution to

$$ (4.4) \qquad \bar{P}^{n+1} + \frac{1}{c_{\lambda}}\bar{P} - \frac{c_rT}{c_{\lambda}} = 0. $$

From now on, we will fix $M > 0$, $c_{\lambda} > 1$, and $n > 0$ and consider the unique stationary solution $\bar{\rho} = \bar{\rho}(c_r)$ (and $\bar{P} = \bar{P}(c_r)$) as a function of the bifurcation parameter $c_r > 0$.

**LEMMA 4.2.** The function $c_r \mapsto \bar{P}(c_r)$ is strictly increasing on $(0, +\infty)$, with $\lim_{c_r \to 0} \bar{P}(c_r) = 0$ and $\lim_{c_r \to +\infty} \bar{P}(c_r) = +\infty$. Moreover, $c_r \mapsto \frac{\bar{P}(c_r)}{c_r}$ is strictly decreasing on $(0, +\infty)$, with $\lim_{c_r \to 0} \frac{\bar{P}(c_r)}{c_r} = T$, and $P(c_r) = o(c_r)$ as $c_r \to +\infty$. Lastly, $0 < \bar{P}'(c_r) < \frac{P(c_r)}{c_r}$ for all $c_r > 0$.

*Proof.* By the Implicit Function Theorem, $\bar{P}$ is differentiable with respect to $c_r$, and by derivating Eq. (4.4) with respect to $c_r$, we get

$$ \bar{P}'(c_r) (c_{\lambda}(n+1)\bar{P}(c_r)^n + 1) = T. $$
Multiplying by $\overline{P}(c_r)$, and using Eq. (4.4), we get

$$ \overline{P}'(c_r) (\mathcal{T}c_r + c_\lambda n \overline{P}(c_r)^{n+1}) = \mathcal{T}\overline{P}(c_r), $$

so

$$ \overline{P}'(c_r) = \overline{P}(c_r) \frac{1}{c_r + \frac{c_\lambda n \overline{P}(c_r)^{n+1}}{\mathcal{T}}} < \frac{\overline{P}(c_r)}{c_r}. $$

Hence $P' > 0$ so $\overline{P}$ is strictly increasing, and by Eq. (4.4), $\lim_{c_r \to +\infty} \overline{P}(c_r) = +\infty$. Moreover,

$$ \left( \frac{\overline{P}(c_r)}{c_r} \right)' = \frac{\overline{P}'(c_r)}{c_r} - \frac{\overline{P}(c_r)}{c_r^2} = \frac{\overline{P}(c_r)}{c_r^2 + \frac{c_\lambda c_r n \overline{P}(c_r)^{n+1}}{\mathcal{T}}} - \frac{\overline{P}(c_r)}{c_r^2} < 0, $$

so the mapping $c_r \mapsto \frac{\overline{P}(c_r)}{c_r}$ is strictly decreasing and positive, hence it admits a limit $l \ge 0$ when $c_r$ goes to $+\infty$. Dividing Eq. (4.4) by $c_r^{n+1}$ and letting $c_r$ go to $+\infty$, we get $l^{n+1} = 0$, so that $l = 0$. Then $\overline{P}(c_r) = o(c_r)$ as $c_r$ goes to $+\infty$. As $\overline{P}$ is positive strictly increasing on $(0, +\infty)$, it admits a nonnegative limit when $c_r$ goes to 0. Letting $c_r$ go to 0 in Eq. (4.4), we get that this limit has to be 0. Furthermore, from Eq. (4.4), we have directly:

$$ \frac{\overline{P}}{c_r} = \frac{\mathcal{T}}{c_\lambda} \frac{1}{\overline{P}^n + \frac{1}{c_\lambda}}, \quad \square $$

$$ \text{so that } \lim_{c_r \to 0} \frac{\overline{P}(c_r)}{c_r} = \mathcal{T}. $$

Now, let us investigate the stability of this stationary solution.

**PROPOSITION 4.3.** We consider the spawning cycle model Eq. (3.1) instantiated with the functions and parameters given by Eq. (4.1)-(4.2). We assume that:

(i) For all $\lambda \in \mathbb{C}$, $c_r > 0$, if $\Delta(\lambda, c_r) = 0$, then $\frac{\partial \Delta(\lambda, c_r)}{\partial \lambda} \neq 0$.

(ii) There exists $c_r^* > 0$ and $\theta^* > 0$ such that $\Delta(i\theta^*, c_r^*) = 0$, $\Re(\frac{\partial \Delta(i\theta^*, c_r^*)}{\partial \lambda}) > 0$, and, for all $\theta \in (0, +\infty) \setminus \{\theta^*\}$, $\Delta(c_r^*, i\theta) \neq 0$.

Then a periodic solution bifurcates from the stationary solution $\bar{\rho}(x) = \frac{\overline{P}(c_r^*)}{\overline{T}f(x)}$, when $c_r = c_r^*$, with $\overline{P}(c_r^*)$ the unique positive solution to Eq. (4.4).

*Proof.* We apply directly Theorem 6.2.7. in [37], which requires to prove transversality conditions 6.2.1.(b), i.e. that there exists a unique pair of conjugated simple roots of $\Delta$ that cross the complex plane imaginary axis as $c_r$ crosses $c_r^*$. First, let us notice that with Eq. (4.4) and condition (i), all roots of $\Delta$ are simple. Furthermore, by Eq. (4.3) and Eq. (4.4):

$$ \begin{aligned} b'(\overline{P}(c_r)) &= -\frac{c_\lambda c_r n \overline{P}(c_r)^{n-1}}{(1 + c_\lambda \overline{P}(c_r)^n)^2} \\ &= -\frac{c_\lambda n}{\overline{T}^2 c_r} \overline{P}(c_r)^{n+1} = -\frac{c_\lambda n}{\overline{T}^2 c_r} \left( \frac{\mathcal{T}c_r}{c_\lambda} - \frac{\overline{P}(c_r)}{c_\lambda} \right) = \frac{n}{\overline{T}^2 c_r} \overline{P}(c_r) - \frac{n}{\overline{T}}. \end{aligned} $$

Then, by Lemma 4.2,

$$ \frac{d}{dc_r} (b'(P)(c_r)) = \frac{n}{\overline{T}^2 c_r} \left( -\frac{\overline{P}(c_r)}{c_r} + \overline{P}'(c_r) \right) < 0. $$
As, for all $\lambda \in \mathbb{C}$, $c_r > 0$, if $\Delta(\lambda, c_r) = 0$, then $\frac{\partial \Delta(\lambda, c_r)}{\partial \lambda} \neq 0$, we can apply the implicit holomorphic function Theorem to $\Delta$ around $(\lambda^* := i\theta^*, c_{r^*})$. This theorem states the existence of $\epsilon > 0$, and a continuously differentiable map $\hat{\lambda}: (c_r^* - \epsilon, c_r^* + \epsilon) \to \mathbb{C}$ such that $\Delta(\lambda, c_r) = 0$ if and only if $\lambda = \hat{\lambda}(c_r)$, for all $(\lambda, c_r) \in \mathcal{B}_{\mathbb{C}}(\lambda^*, \epsilon) \times (c_r^* - \epsilon, c_r^* + \epsilon)$. Moreover, as $\Delta(\hat{\lambda}(c_r), c_r) \equiv 0$ for all $c_r \in (c_r^* - \epsilon, c_r^* + \epsilon)$, we have by Lemma 4.5:

$$ \frac{d\hat{\lambda}(c_r)}{dc_r} = - \frac{1}{\frac{\partial \Delta(\hat{\lambda}(c_r), c_r)}{\partial \lambda}} \frac{\partial \Delta(\hat{\lambda}(c_r), c_r)}{\partial c_r} = \frac{1}{\frac{\partial \Delta(\hat{\lambda}(c_r), c_r)}{\partial \lambda}} \frac{db'(\bar{P})(c_r)}{dc_r} \frac{1}{b'(\bar{P})(c_r)}. $$

In particular, at $c_r = c_{r^*}$, we obtain

$$ \frac{d\hat{\lambda}(c_r^*)}{dc_r} = \frac{1}{\frac{\partial \Delta(\hat{\lambda}(c_r^*), c_r^*)}{\partial \lambda}} \frac{1}{b'(\bar{P})(c_r^*)} \frac{db'(\bar{P})(c_r^*)}{dc_r}, $$

so,

$$ \frac{d\Re(\hat{\lambda}(c_{r^*}))}{dc_r} > 0. \quad \square $$

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

Now, to verify explicitly the assumptions of Proposition 4.3, we need to specify function $f$. We will assume that $f \equiv 1$ in order to simplify the calculations. Explicit calculations become indeed very tedious even with simple choices such as an affine function $f$. With some direct calculations detailed in Appendix E, we can prove the following proposition.

**PROPOSITION 4.4.** If $f \equiv 1$, the characteristic function $\Delta$ defined by Eq. (3.13) can be expressed as:

$$ \Delta(\lambda, c_r) = 1 - \frac{2b'(\bar{P})(c_r)M}{\lambda^2} (e^{-\frac{\lambda}{2}} - 1)^2, \quad \forall \lambda \in \mathbb{C} \setminus \{0\} $$

$$ \text{and } \Delta(0, c_r) = 1 - \frac{1}{2} b'(\bar{P})(c_r) M > 0. $$

Then $\lambda$ is a root of $\Delta(\cdot, c_r)$ if and only if $\lambda \neq 0$ and

$$ (4.5) \qquad \frac{e^{-\frac{\lambda}{2}} - 1}{\lambda} = \pm i \sqrt{\frac{-1}{2b'(\bar{P})(c_r)M}}. $$

Let $k \in \mathbb{N}$ and $c_{r,k} > 0$ verifying

$$ (4.6) \qquad k\pi + \frac{\pi}{2} = \sqrt{n}\sqrt{1 - \frac{2}{M}\frac{\bar{P}(c_{r,k})}{c_{r,k}}}, $$

then $\Delta(\cdot, c_{r,k})$ has a unique pair of conjugated purely imaginary roots $(\lambda_k, \bar{\lambda}_k)$ with $\lambda_k = (4k\pi + 2\pi) i$. For all $c_r \neq c_{r,k}$, $\Delta(\cdot, c_r)$ does not admit any purely imaginary root.

Let us remark that by Lemma 4.2, $c_r \mapsto 2\sqrt{n}\sqrt{1 - \frac{2}{M}\frac{\bar{P}(c_r)}{c_r}}$ is a strictly increasing continuous function on $(0, +\infty)$, with $\lim_{c_r \to 0} 2\sqrt{n}\sqrt{1 - \frac{2}{M}\frac{\bar{P}(c_r)}{c_r}} = 0$, and

$$ \lim_{c_r \to +\infty} 2\sqrt{n}\sqrt{1 - \frac{2}{M}\frac{\bar{P}(c_r)}{c_r}} = 2\sqrt{n}. \text{ So, by Proposition 4.5, } \Delta \text{ admits } K + 1 \text{ purely} $$
imaginary roots $(\lambda_k)_{0 \le k \le K}$ with $K = \lfloor \frac{\sqrt{n}}{\pi} - \frac{1}{2} \rfloor$. In particular, for $n \le \frac{\pi^2}{4}$, $\Delta$ does not admit any purely imaginary root, whatever the value of $c_r$, and, for $n > \frac{\pi^2}{4}$, it exists at least one $c_r$ such that $\Delta(\cdot, c_r)$ admits a pair of purely imaginary roots. By propositions 3.5 and 4.4, we can easily deduce the next Lemma.

LEMMA 4.5. The derivative of $\Delta$ with respect to $\lambda$ is given by

$$ \frac{\partial \Delta(\lambda, c_r)}{\partial \lambda} = \frac{2b'(\bar{P})(c_r)M}{\lambda^2} \left( e^{-\lambda} - e^{-\frac{\lambda}{2}} \right) + \frac{2}{\lambda} (1 - \Delta(\lambda)). $$

If $\lambda$ is a root of $\Delta$, then

$$ \frac{\partial \Delta(\lambda, c_r)}{\partial \lambda} = \frac{2b'(\bar{P})(c_r)M}{\lambda^2} \left( e^{-\lambda} - e^{-\frac{\lambda}{2}} \right) + \frac{2}{\lambda} = \frac{1}{1 - e^{\frac{\lambda}{2}}} + \frac{2}{\lambda} \neq 0. $$

In particular, all roots of $\Delta$ are simple. If $\lambda$ is a purely imaginary root, then

$$ \frac{\partial \Delta(\lambda)}{\partial \lambda} = \frac{1}{2} - \frac{i}{\sqrt{-2b'(\bar{P})(c_r)M}}. $$

Moreover, the derivative of $\Delta$ with respect to $c_r$ at $(c_r, \lambda)$ where $\lambda$ is a root of $\Delta$ is given by

$$ \frac{d\Delta(\lambda, c_r)}{dc_r} = - \frac{db'(\bar{P})(c_r)}{dc_r} \frac{1}{b'(\bar{P}(c_r))}. $$

Therefore, we can apply Proposition 4.3 to state the main result of this subsection.

PROPOSITION 4.6. Let us consider the spawning cycle model Eq. (3.1) with model parameters given by Eq. (4.1)-(4.2), and $f \equiv 1$. For any $M > 0$, $c_\lambda > 1$, and $n > \frac{\pi^2}{4}$, there is a Hopf bifurcation at $c_r = c_{r,k}$, with $k \in \{0, \dots, \lfloor \frac{\sqrt{n}}{\pi} - \frac{1}{2} \rfloor\}$, where $c_{r,k}$ verifies Eq. (4.6). In particular, a periodic solution bifurcates from the stationary solution $\bar{\rho}(c_{r,k}) = \frac{2\bar{P}(c_{r,k})}{M}$ when $c_r = c_{r,k}$, with $\bar{P}(c_{r,k})$ the unique positive solution to Eq. (4.4).

### 4.1.2. Numerical simulations.
We now provide a numerical illustration of the Hopf bifurcation Proposition 4.6 through direct simulations of the spawning cycle model (see Appendix SM1 for details on the numerical scheme). We set $f \equiv 1$, $M = 1$, $n = 4$, $c_\lambda = 2$. From Proposition 4.6, we have a unique Hopf bifurcation at $c_r^0 \approx 4.944$, obtained by solving numerically Eq. (4.6) (see the bifurcation diagram on Figure 3). On the top left panel of Figure 2, we have set $c_r = 4 < c_r^0$. The unique positive solution to Eq. (4.4) is $\bar{P} \approx 0.889$, and the associated stationary solution is given by $\bar{\rho}(x) = \frac{2\bar{P}}{M} \approx 1.778$. Starting from the initial condition $\rho_0(x) \equiv 0$, we observe that $P$ converges to $\bar{P}$, so that the population of large-sized oocytes converges to the associated stationary solution $\bar{\rho}$, as expected. In the remaining panels, we have set $c_r = 10 > c_r^0$. Still starting from $\rho_0(x) \equiv 0$, we observe that the population of large-sized oocytes converges to a periodic solution, which seems to be stable. Indeed, on the top right panel ($P$ with respect to time) and bottom left panel (heatmap of the density $\rho$), we can see the persistence of a periodic pattern. The shape of the density is shown on the bottom right panel. We observe a single-peaked density that propagates rightwards with time. It is worth noting that the shape of this periodic
solution is similar to the large-sized oocyte size distribution observed in experimental data (see, e.g., Figure 16 in [46]).

<figure><img src="image_4.png" /><figcaption>FIG. 2. Illustration of the Hopf bifurcation and a periodic solution to the spawning cycle model, with $f \equiv 1$, $n = 4$, $M = 1$, $c_{\lambda} = 2$. Top left panel shows the dynamics of P with $c_r = 4$ while top right panel shows the dynamics of P with $c_r = 10$. Bottom panels show the dynamics of the periodic solution $\rho$ for $c_r = 10$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3. Hopf bifurcation diagram with respect to $c_r$, for the spawning cycle model with $f \equiv 1$, $n = 4$, $M = 1$, $c_{\lambda} = 2$. The stationary solution branch is computed by solving Eq. (4.4) for several $c_r$ on a grid of $c_r$ values, and the branches corresponding to the amplitude of the periodic solution are determined numerically from simulated solutions.</figcaption></figure>

We keep $f \equiv 1$, $M = 1$, $n = 4$, $c_{\lambda} = 2$, but we now perform a numerical simulation with a non monotonous time-dependent basal rate of newly recruited large-sized oocytes $c_r(t)$, such that:

$$c_r(t) = \begin{cases} 10(1 - e^{-0.1t}) & \text{if } t \le 60 \\ 10(1 - e^{-6})e^{-0.1(t-60)} & \text{if } t > 60. \end{cases}$$

The expression of $c_r(t)$ roughly mimics the change in the basal recruitment rate along aging (see left panel of Figure 4). The corresponding P dynamics (which fully characterizes $\rho$ dynamics) are plotted on the right panel of Figure 4. One can observe a 3 step pattern. For small $t$, hence in early (prepubertal) life, $c_r$ is lower than the bifurcation value $c_r^0$, there is a single stable stationary solution. Then, during the mid-life (reproductive) period, the cyclic solution appears and is maintained till the value of
$c_r$ drops below the critical value (senescence). These results seem to indicate that there may be a threshold on the number of recruited oocytes to trigger and maintain the spawning cycle. This could explain why some elder individuals are no longer able to spawn.

Note that we observe a delay to the bifurcation, which occurs as the current point lies in the vicinity of the unstable stationary solution, as classically encountered in slow-fast dynamical systems.

<figure><img src="image_6.png" /><figcaption>FIG. 4. Simulation of the spawning cycle model with $n = 4$, $M = 1$, $c_λ = 2$ (right panel), with a time-dependent basal recruitment $c_r$ (left panel).</figcaption></figure>

## 4.2. A model of small-sized oocyte dynamics (previtello).
We now focus on the population of small-sized (previtellogenic) oocytes, whose slower dynamics drive the whole oocyte population dynamics on the lifespan. For this application, the lower size boundary ($x = 0$) in Eq. (3.1) now corresponds to the size of newly recruited growing oocytes, and the greater size boundary $x = 1$ to the size of oocytes entering the vitellogenic stage. In the previtellogenic phase, AMH both slows down oocyte growth and exerts a critical inhibitory effect on oocyte recruitment, which appears to be associated with an upstream inhibition of germline cell proliferation [48, 49].

Accordingly, we formulate the recruitment and growth rates as follows:

$$ (4.7) \qquad \Lambda(P, x) = f(x)g(P) = \frac{f(x)}{1+c_{\lambda}P^n}, \quad r(P) = \frac{c_r}{1+P^n}, \quad n > 1, $$

so that $b = \frac{r}{g} = \frac{c_r(1+c_\lambda P^n)}{1+P^n}$. The parameter $c_r > 0$ is the basal rate of newly recruited growing (previtellogenic) oocytes, and parameter $c_\lambda > 0$ tunes the intensity of the feedback exerted onto the growth rate. The function $f$ is strictly positive continuously differentiable and $\omega$ is non negative Lipschitz continuous. Like for the spawning cycle model, we will study the asymptotic behavior of $\rho$ with respect to $c_r$.

Let us note that this model has a similar asymptotic behavior as that with $g(P) = 1 + P^n$, and $r(P) = c_r(1 + c_\lambda P^n)$, where there is a double positive feedback on the growth rate and recruitment rate. Such a formulation would be well-suited for representing the estradiol-controlled dynamics of oocytes in the vitellogenic phase (see Figure 1).

### 4.2.1. Bifurcation analysis.
By application of Proposition 3.1, we know that the stationary solutions are given by $\bar{\rho}(c_r)(x) = \frac{b(\bar{P})(c_r)}{f(x)}$ where $\bar{P}(c_r)$ is a (not necessarily unique) positive solution to

$$ (4.8) \qquad \bar{P}(c_r) = \mathcal{T}c_r \frac{1+c_\lambda \bar{P}(c_r)^n}{1+\bar{P}(c_r)^n}, $$
with $\mathcal{T} = \int_0^1 \frac{\omega(x)}{f(x)}dx$. We enunciate the stability properties of these stationary solutions in the next proposition, which is illustrated in Figures 5 and 6.

**PROPOSITION 4.7.** For each $c_r > 0$, $c_\lambda > 0$, $n > 1$, $\omega$ a nonnegative Lipschitz continuous function, f a strictly positive continuously differentiable function, the small-sized oocyte model, Eq. (3.1) with model parameters given by Eq. (4.7), has the following stationary solutions:

* If $c_\lambda < 1$, then for each $c_r > 0$ there is a unique stationary solution.

* If $c_\lambda \ge 1$ and the following condition is satisfied:

$$ (4.9) \qquad C := ((1-n)c_\lambda + 1 + n)^2 - 4c_\lambda > 0, $$

we define

$$ P_{\pm} = \left( \frac{(n-1)c_{\lambda} - 1 - n \pm \sqrt{C}}{2c_{\lambda}} \right)^{\frac{1}{n}}, $$

and

$$ (4.10) \qquad c_{r\pm} = \frac{P_{\mp}}{\mathcal{T}} \frac{1 + P_{\mp}^n}{1 + c_{\lambda} P_{\mp}^n}. $$

Then, we have the following:

- if $c_r < c_{r-}$ there is a unique locally asymptotically stable stationary solution $\overline{\rho}_1(c_r)$,

- if $c_r = c_{r-}$ there are two stationary solutions $\overline{\rho}_1(c_r) < \overline{\rho}_{23}(c_r)$. $\overline{\rho}_1(c_r)$ is locally asymptotically stable, and $\overline{\rho}_3(c_r)$ is non hyperbolic,

- if $c_{r-} < c_r < c_{r+}$, then there are three stationary solutions $\overline{\rho}_1(c_r) < \overline{\rho}_2(c_r) < \overline{\rho}_3(c_r)$. $\overline{\rho}_1(c_r)$ and $\overline{\rho}_3(c_r)$ are locally asymptotically stable, and $\overline{\rho}_2(c_r)$ is unstable,

- if $c_r = c_{r+}$, there are two stationary solutions $\overline{\rho}_{12}(c_r) < \overline{\rho}_3(c_r)$. $\overline{\rho}_{12}(c_r)$ is non hyperbolic, and $\overline{\rho}_3(c_r)$ is locally asymptotically stable,

- if $c_r > c_{r+}$, then there is a unique locally asymptotically stable stationary solution $\overline{\rho}_3(c_r)$.

To sum up, when $c_\lambda \ge 1$ is fixed and Eq. (4.9) is satisfied, two successive saddle-node bifurcations occur at $c_r = c_{r-}$ and $c_r = c_{r+}$.

* If $c_\lambda \ge 1$ and condition (4.9) is not satisfied, then for each $c_r$ there is a unique locally asymptotically stable stationary solution.

Condition Eq. (4.9) implies that $b$ is an increasing Hill function with a sufficient steep slope (see Figure 5).

*Proof.* If $c_\lambda < 1$, then $b$ is strictly decreasing, so that Eq. (4.8) admits a unique positive solution. Else, we can easily check that Eq. (4.8) can admit 1,2 or 3 positive solutions, depending on the shape of $b$ (see Figure 5). If Eq. (4.9) is satisfied, then the critical values of $c_r$ (transition from 1 to 3 stationary solutions), for which we have two stationary solutions, are given by the following equations:

$$ (4.11) \qquad \overline{P} = \mathcal{T}c_r \frac{1 + c_\lambda \overline{P}^n}{1 + \overline{P}^n} \quad \text{and} \quad 1 = \mathcal{T}b'(\overline{P}) = \mathcal{T}c_r \frac{n \overline{P}^{n-1} (c_\lambda - 1)}{\left(1 + c_\lambda \overline{P}^n\right)^2}. $$
Positive solutions $\bar{P}_{\pm}$ to Eq. (4.11) can be determined with few calculations, as in [35], and used to deduce $c_{r\pm}$ given by Eq. (4.10). Finally, when $c_{\lambda} \ge 1$ and Eq. (4.9) is not satisfied, Eq. (4.11) has no solution, hence Eq. (4.8) admits always one unique positive solution.

Now, concerning stability results, when $c_{\lambda} \ge 1$, $b$ is an increasing function with respect to $P$, so that we can use Corollary 3.7 to deduce the local stability of stationary solutions. More precisely, Corollary 3.7 implies that $\bar{\rho}$ (associated with $\bar{P}$) is locally asymptotically stable if the slope of the tangent to the sigmoid $P \mapsto Tb(P)$ at point $\bar{P}$ is strictly greater than the slope of $P \mapsto P/c_r$ (which is $1/c_r$). Conversely, if the slope of the tangent curve is strictly smaller than $1/c_r$, then $\bar{\rho}$ is unstable. Therefore, we can easily deduce the stability of stationary solutions (see Figure 5 for an illustration of the proof). □

It should be noted that when $c_{\lambda} < 1$, then $b$ is decreasing, leading possibly to qualitatively similar phenomena as in the spawning cycle application (subsection 4.1).

<figure><img src="image_7.png" /><figcaption>FIG. 5. Illustration of Proposition 4.7, with $n = 6$, $f(x) = 1 + 2x^2$, $\omega(x) = x$. Dotted lines are graphical representations of $P \mapsto \frac{P}{c_r}$, black curves of $P \mapsto Tb(P)$, and black points are intersection points. In A, B, C panels, $c_{\lambda} = 5$, such that we are in the case $c_{\lambda} \ge 1$ and condition Eq. (4.9) holds true. In A $c_r < c_{r-}$, in B $c_{r-} < c_r < c_{r+}$, and in C $c_r > c_{r+}$. In panel D, $c_{\lambda} = 1.8$, such that condition Eq. (4.9) is not satisfied.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. The left panel shows the bifurcation diagram with respect to $c_λ$ and $c_r$ of the small-sized oocyte model Eq. (4.7) with $n = 6$, $f(x) = 1+2x^2$, $\omega(x) = x$. Captioned areas A, B, C, D correspond to the cases illustrated in Figure 5. The right panel represents the bifurcation diagram with respect to $c_r$ of the small-sized oocyte model Eq. (4.7) with $n = 6$, $f(x) = 1+2x^2$, $\omega(x) = x$, and with $c_λ$ fixed to 5.</figcaption></figure>

### 4.2.2. Numerical simulation with a time-dependent recruitment rate.

Let us consider the small-sized oocyte model with $n = 6$, $f(x) = 1+2x^2$, $\omega(x) = x$, and $c_λ = 5$. The choice of $\omega$ is based (roughly) on literature about AMH expression [39]. AMH expression increases during the previtellogenic phase until reaching a peak at the late previtellogenic stage, hence the parsimonious choice of a linear function for $\omega$. $f$ is an increasing quadratic function with respect to $x$, accounting for the acceleration of the intrinsic grow rate in this size range. $n$ is fixed arbitrarily and $c_λ$ as well, but in such a way that $b$ is increasing. The basal recruitment rate $c_r$ is slowly increasing, and then slowly decreasing with respect to time (top left panel of Figure 7), as it can be expected at the scale of a whole fish life:

$$c_r(t) = \begin{cases} 0.75 + 1.5 \sin \left( \frac{\pi t}{10000} \right) & \text{if } 0.75 + 1.5 \sin \left( \frac{\pi t}{10000} \right) < 2, \\ 2 & \text{if } 0.75 + 1.5 \sin \left( \frac{\pi t}{10000} \right) \ge 2. \end{cases}$$

We observe an hysteresis loop (bottom panel of Figure 7), with first the transition from a stationary solution branch with few oocytes to a stationary solution branch with significantly more oocytes when $c_r$ increases, and then the reverse transition when $c_r$ decreases. We recall that the stationary solutions have all the same shape (proportional to $1/f(x)$), so that the transition from one branch to another simply corresponds to a change in the size of the population, not to a change in the shape of the density.
<figure><img src="image_9.png" /><figcaption>FIG. 7. Simulation of the small-sized oocyte model with $n = 6$, $f(x) = 1 + 2x^2$, $\omega(x) = x$, $c_\lambda = 5$ and a time-dependent basal recruitment rate $c_r$ as shown in the top left panel. The top right panel shows the corresponding dynamics of $P$ with respect to $t$. On the bottom panel, the trajectory $P(t)$ (orange curve) is overlaid on the bifurcation diagram (black curve).</figcaption></figure>

**5. Conclusion.** In this manuscript, we have introduced a fairly general quasi-linear size-structured population model defined by Eq. (1.1)-(1.2) and established the existence and uniqueness of global solutions in Theorem 2.2. The originality of the model mainly comes from the inflow boundary condition, which is formulated as a recruitment term rather than as a birth term. Our main mathematical contribution lies in the long-time analysis of the model under some specific assumptions. As a general theory for generic quasilinear size-structured population models seems out of reach, we have limited our study to model Eq. (3.1)-(1.2), assuming that the growth rate is separable (the feedback onto the growth of individuals is homogeneous across the population), and that there is no death. We have shown that this model is asymptotically equivalent to the semilinear model defined by Eq. (3.4).

The long-time analysis has revealed the richness of the model, which exhibits monostability, multistability (see Corollary 3.9 and Proposition 4.7) and oscillatory solutions (see Proposition 4.6), depending on the recruitment-to-growth ratio. When the ratio is a decreasing function of the population weighted average (which typically indicates the presence of negative feedback onto recruitment and/or positive feedback on growth rate), we have shown the existence of a Hopf bifurcation, and confirmed numerically the stability of the oscillatory solution. When the ratio is an increasing sigmoidal function of the population weighted average, we have demonstrated the existence of two successive (supercritical and subcritical) saddle-node bifurcations leading to possible bistability. However, as we considered only one interaction term (i.e. $N = 1$) in Sec. 4, we did not investigate the case of a recruitment-to-growth ratio involving different directions of variation in the different interaction terms. It could be the matter of future work.

From an application perspective, this work represents a first step towards understanding the mechanisms controlling oocyte dynamics in fish, particularly the dynamics of small-sized oocyte recruitment, which remains only partially understood.
In future works, it would be interesting to embed the submodels representing separately the main steps of oogenesis into a comprehensive model encompassing all stages of oocyte growth, and to study the influence of the recruitment on the whole population dynamics.

The assumption of a separable $f(x)g(P)$ growth rate is well suited to the submodels, since the scalar $P(t)$ term in each submodel can be interpreted as the level of the prominent hormone, which both emanate from, and act similarly on a specific and rather homogenous oocyte subpopulation. More generally, the separable formulation corresponds to situations in which the nonlocal control – hormonal control in our case– acts in a uniform way according to the structuring variable – size as a marker of maturity in our case –. In some sense, the $g(P(t))$ term amounts to a time-varying scaling term modulating the intrinsic part of the growth rate driven by the $f(x)$ term.

The separable formulation would not be as valid for the complete model spanning the whole oogenesis, involving altogether AMH, Estradiol and DHP in a vector $P(t)$, since for instance immature stages are rather insensitive to the Estradiol-mediated and even more DHP-mediated control terms. Generally speaking, the separability assumption would fail in the case when parts of the population would be insensitive to the nonlocal control, or if this control could have antagonistic effects depending on the size/maturity level.

Coupling the submodels can be achieved either through the flux boundary conditions, or through the design of a single PDE with a velocity that depends both on $x$ and $P$ (quasilinear PDE). To study the behavior of such a model, numerical methods such as the pseudospectral approximation coupled with numerical tools for bifurcation analysis of ordinary differential equations would be useful [7, 40].

## Appendix A. (Local) Existence and Uniqueness.

It is easy to see that Eq. (1.1) with a given $P \in C([0, T], (\mathbb{R}^+)^N)$ admits a unique (weak) solution given by the following explicit formula on $[0, T] \times [0, 1]$:

$$ egin{cases} 
\rho(t, x) = \rho_0(\beta_P(t, x))\mathcal{F}_P(t; 0, \beta_P(t, x))\frac{\partial\beta_P}{\partial x}(t, x), & (t, x) \in D_\beta \\ 
\rho(t, x) = -r(\alpha_P(t, x), P(\alpha_P(t, x)))\mathcal{F}_P(t; \alpha_P(t, x), 0)\frac{\partial\alpha_P}{\partial x}(t, x), & (t, x) \in D_\alpha. 
\end{cases} $$

This standard result can be obtained by the method of characteristics. Next Propositions A.1-A.4 will ensure that the corresponding result holds true in the nonlinear case (see Theorem 2.2), i.e. for $P(t) = \int_0^1 \omega(x)\rho(t, x)dx$.

<figure><img src="image_10.png" /><figcaption>FIG. 8. The characteristic $X_{0,P}$ starting from the origin split the domain in two parts: $D_\alpha$ (in salmon) and $D_\beta$ (in blue). Characteristics in $D_\alpha$ intersect the time axis at time $\alpha$, while characteristics in $D_\beta$ intersect the size axis at size $\beta$.</figcaption></figure>
To ease the formulation of stability estimates, we introduce the following notations, where $T > 0$ is a fixed time:

$$ (A.1) \qquad M := \max \left( 1 + \|\rho_0\|_{\infty} \|\omega\|_{\infty}, N \|\omega\|_{\infty} (\|\rho_0\|_{\infty} + T \|r_0\|_{\infty}) \exp \left( TN \|\omega\|_{\infty} \sup_{1 \le i \le N} \|r_i\|_{\infty} \right) \right), $$

and

$$ (A.2) \qquad \delta := \min \left( T, \frac{1}{\|\omega\|_{\infty} \left( \|r_0\|_{\infty, [0,T]} + M \sum_i \|r_i\|_{\infty, [0,T]} \right)}, \frac{1}{2L_{\Lambda, [0,M]^N \times [0,1]}}, \frac{1}{c_1 + c_2 \|\rho_0\|_{\infty}} \right). $$

where constants $c_1, c_2$, depend on $T$, $\|r\|_{\infty, [0,T] \times [0,M]^N}$, $\|\Lambda\|_{\infty, [0,M]^N \times [0,1]}$, and Lipschitz constants $L_{r, [0,M]^N \times [0,1]}$, $L_{\Lambda, [0,M]^N \times [0,1]}$, $L_{\mu, [0,M]^N \times [0,1]}$, $L_{\omega}$, and can be explicitly deduced from Eq. (A.5).

**PROPOSITION A.1.** Let $T > 0$, and $\rho$ a weak solution to Eq. (1.1) up to time $T$. Therefore, for all $\tau \in [0, T]$ and $\varphi \in C^1([0, \tau] \times [0, 1])$ such that $\varphi(t, 1) = 0$ for all $t \in [0, \tau]$, we have:

$$ \int_0^\tau \int_0^1 \rho(t, x) [\partial_t \varphi(t, x) + \Lambda(P(t), x) \partial_x \varphi(t, x) - \mu(P(t), x) \varphi(t, x)] dxdt \\ + \int_0^\tau r(t, P(t)) \varphi(t, 0) dt + \int_0^1 \rho_0(x) \varphi(0, x) dx - \int_0^1 \rho(\tau, x) \varphi(\tau, x) dx = 0. $$

*Proof.* The proof of this proposition is similar to that provided for Lemma 2.2. in [11]. □

**PROPOSITION A.2.** Let $\rho$ be a weak solution to Eq. (1.1) up to $\delta$, with $\delta$ verifying (A.2). Then $\rho$ verifies the quasi-explicit formulation Eq. (2.5).

*Proof.* We consider $0 \le \tau \le \delta$, and $\psi$ the solution to the following (linear) dual problem:

$$ \left\{ \begin{array}{l} \partial_t \psi(t, x) + \Lambda(P(t), x) \partial_x \psi(t, x) - \mu(P(t), x) \psi(t, x) = 0, \quad x \in [0, 1], \quad t \in [0, \tau] \\ \psi(t, 1) = 0, \quad t \in [0, \tau] \\ \psi(\tau, x) = \psi_0(x), \quad x \in [0, 1] \end{array} \right. $$

with $\psi_0 \in C_c^1(0, 1)$. By the characteristic method, $\psi$ verifies:

$$ (A.3) \qquad \left\{ \begin{array}{l} \psi(t, x) = \psi_0(X_P(\tau; t, x)) \mathcal{F}_P(\tau; t, x), \quad 0 \le x \le X_P(t, \tau, 1), \quad 0 \le t \le \tau \\ \psi(t, x) = 0, \quad X_P(t, \tau, 1) \le x \le 1, \quad 0 \le t \le \tau. \end{array} \right. $$

As $\psi_0 \in C_c^1(0, 1)$, we have $\psi_0(1) = \psi'_0(1) = 0$, and thus $\psi \in C^1([0, \tau] \times [0, 1])$. Then by applying Proposition A.1 with $\varphi = \psi$, we get:

$$ \int_0^1 \rho(\tau, x) \psi_0(x) dx = \int_0^1 \rho_0(x) \psi(0, x) dx + \int_0^\tau r(t, P(t)) \psi(t, 0) dt. $$
With (A.3) and by doing the change of variables $y = X_P(\tau; 0, x)$ (see Proposition 2.1), the first r.h.s. (right-hand side term) becomes:

$$ \int_0^1 \rho_0(x)\psi(0,x)dx = \int_{X_{0,P}(\tau)}^1 \rho_0(\beta_P(\tau, y))\mathcal{F}_P(\tau; 0, \beta_P(\tau, y)) \frac{\partial \beta_P}{\partial x}(\tau, y)\psi_0(y)dy. $$

Additionally, with (A.3) and the change of variables $x = X_P(\tau; t, 0)$ (see Proposition 2.1), the second r.h.s. becomes:

$$ \int_0^\tau r(t, P(t))\psi(t, 0)dt = \int_0^{X_{0,P}(\tau)} -r(\alpha_P(\tau, x), P(\alpha_P(\tau, x)))\mathcal{F}_P(\tau; \alpha_P(t, x), 0) \frac{\partial \alpha_P}{\partial x}(\tau, x)\psi_0(x)dx. $$

Since these equalities are verified for every $\psi_0 \in C_c^1(0, 1)$ and $0 \le \tau \le \delta$, then $\rho$ verifies the quasi-explicit formulation Eq. (2.5) on $[0, 1] \times [0, \delta]$. $\square$

The following result is a corollary of Proposition A.2.

**PROPOSITION A.3 (Local uniqueness).** Let $\delta > 0$ verifying (A.2), then Eq. (1.1) admits at most one weak solution up to time $\delta$.

*Proof.* Let us suppose that there exist two weak solutions $\rho$ and $\tilde{\rho}$. Thanks to Proposition A.2, we know that $\rho$ verifies the quasi-explicit formulation Eq. (2.5) with $P = \int_0^1 \omega(x)\rho(t, x)dx$, and $\tilde{\rho}$ verifies the quasi-explicit formulation Eq. (2.5) with $P$ replaced by $\tilde{P} = \int_0^1 \omega(x)\tilde{\rho}(t, x)dx$. $P$ and $\tilde{P}$ are both fixed points of the mapping $F$ defined by Eq. (2.7). By Lemma 2.4, we know that $F$ admits a unique fixed point on $\Omega_{\delta,M} = \{P \in C^0([0, \delta])^N : \forall 1 \le i \le N, P_i \ge 0 \text{ and } \|P_i\|_{C^0([0, \delta])} = \sup_{t \in [0, \delta]} P_i(t) \le M\}$, with $M$ defined by (A.1). With the same calculations that led to Eq. (2.9), we can state (by keeping the same notations of constants $B$ and $C$), that

$$ \sup_{t \in [0, \delta]} P_i(t) \le Be^{C\delta} \le Be^{CT} \le M \quad \text{and} \quad \sup_{t \in [0, \delta]} \tilde{P}_i(t) \le Be^{C\delta} \le Be^{CT} \le M. $$

Hence $P$ and $\tilde{P}$ belong to $\Omega_{\delta,M}$, and by uniqueness of the fixed point $P \equiv \tilde{P}$ on $[0, \delta]$, and thus $\rho = \tilde{\rho}$.

**PROPOSITION A.4 (Local existence).** Let $0 < \delta < T$ with $\delta$ verifying (A.2), $M > 0$ verifying (A.1) and $P$ a fixed point of $F$ defined by Eq. (2.7). Then $\rho$ defined by Eq. (2.5) up to time $\delta$ is a weak solution to Eq. (1.1).

*Proof.* Firstly, let us verify that $\rho$ defined by Eq. (2.5) up to time $\delta > 0$ is in $C^0([0, \delta]; L^1(0, 1)) \cap L^\infty((0, \delta) \times (0, 1))$.

With Eq. (2.5) and regularity assumptions (H2)-(H6) combined with the fact that $P \in \Omega_{\delta,M}$, we obtain the following upper bound:

$$ \sup_{(t,x) \in [0,1] \times [0,\delta]} |\rho(t,x)| \le e^{T\|\Lambda_x\|_{\infty, [0,M]^N} \times [0,1]} \max \left\{ \|\rho_0\|_{\infty}, \frac{\|r\|_{\infty, [0,T] \times [0,M]^N}}{\inf_{[0,M]^N \times [0,1]} \Lambda} \right\}. $$

It remains to show that $\rho \in C^0([0, \delta]; L^1(0, 1))$. Let $0 \le t \le \tilde{t} \le \delta$, and let us prove that $\lim_{\tilde{t} \to t^+} \|\rho(\tilde{t}, \cdot) - \rho(t, \cdot)\|_{L^1(0,1)} = 0$. The case $\tilde{t} \le t$ works in the same way.

For sake of legibility, we will note $X_P(\theta) := X_P(\theta; t, x)$, $\tilde{X}_P(\theta) := X_P(\theta; \tilde{t}, x)$, and
$\alpha_P := \alpha_P(t, x)$, $\tilde{\alpha}_P := \alpha_P(\tilde{t}, x)$, $\beta_P := \beta_P(t; x)$, $\tilde{\beta}_P := \beta_P(\tilde{t}; x)$.
Let $x \in [0, X_{0,P}(t)]$, then with Eq. (2.5) and as $t \le \tilde{t}$:

$$ \rho(\tilde{t}, x) = -r(\tilde{\alpha}_P, P(\tilde{\alpha}_P))\mathcal{F}_P(\tilde{t}; \tilde{\alpha}_P, 0) \frac{\partial \alpha_P}{\partial x}(\tilde{t}, x) \\ = r(\tilde{\alpha}_P, P(\tilde{\alpha}_P))e^{-\int_{\tilde{\alpha}_P}^{\tilde{t}} \mu(P(\theta), X_P(\theta, \tilde{\alpha}_P, 0)) d\theta} \frac{\partial \alpha_P}{\partial x}(\tilde{t}, x). $$

We recall that $P$ is a fixed point of $F$ and in particular $P \in \Omega_{\delta,M}$. Additionally, Definition-Proposition 2.1 gives regularity properties on $\alpha_P$, ensuring that $\tilde{\alpha}_P \xrightarrow{\tilde{t} \to t^+} \alpha_P$ and $\partial_x \alpha_P(\tilde{t}, x) \xrightarrow{\tilde{t} \to t^+} \partial_x \alpha_P(t, x)$. By combining these results with regularity assumptions (H2)-(H6), one can easily check that $\rho(\tilde{t}, \cdot)$ converges pointwise to $\rho(t, \cdot)$ on $[0, X_{0,P}(t)]$. Thanks to Eq. (A.4), we can apply the dominated convergence Theorem which leads to:

$$ \lim_{\tilde{t} \to t^+} \|\rho(\tilde{t}, \cdot) - \rho(t, \cdot)\|_{L^1(0, X_{0,P}(t))} = 0. $$

Furthermore, still with Eq. (2.5)

$$ \int_{X_{0,P}(\tilde{t})}^{1} |\rho(\tilde{t}, x) - \rho(t, x)| dx \le \int_{X_{0,P}(\tilde{t})}^{1} |\rho_0(\tilde{\beta}_P) - \rho_0(\beta_P)| \mathcal{F}_P(\tilde{t}, 0, \tilde{\beta}_P) |\partial_x \beta_P(\tilde{t}, x)| dx \\ + \int_{X_{0,P}(\tilde{t})}^{1} \rho_0(\beta_P) |\mathcal{F}_P(t, 0, \beta_P) \partial_x \beta_P(t, x) - \mathcal{F}_P(\tilde{t}, 0, \tilde{\beta}_P) \partial_x \beta_P(\tilde{t}, x)| . $$

Like before, with Definition-Proposition 2.1, regularity assumptions (H2)-(H6) and the dominated convergence Theorem, one can check that the second right-hand side term goes to 0 when $\tilde{t}$ goes to $t$. It remains to show that the first right-hand side term also vanishes. As $C^1([0, 1])$ is dense in $L^1(0, 1)$, there exists $(\rho_0^n)_{n \ge 0} \in C^1([0, 1])^N$ such that $\rho_0^n \xrightarrow{n \to +\infty} \rho_0$. By doing an upper bounding analogously to Eq. (A.4) and denoting $C$ a constant that does not depend on neither $t$ nor $x$, we have

$$ \int_{X_{0,P}(\tilde{t})}^{1} |\rho_0(\tilde{\beta}_P) - \rho_0(\beta_P)| \mathcal{F}_P(\tilde{t}, 0, \tilde{\beta}_P) |\partial_x \beta_P(\tilde{t}, x)| dx \\ \le C \int_{X_{0,P}(\tilde{t})}^{1} |\rho_0(\tilde{\beta}_P) - \rho_0(\beta_P)| dx, $$

and

$$ \int_{X_{0,P}(\tilde{t})}^{1} |\rho_0(\tilde{\beta}_P) - \rho_0(\beta_P)| dx \le \int_{X_{0,P}(\tilde{t})}^{1} |\rho_0(\tilde{\beta}_P) - \rho_0^n(\tilde{\beta}_P)| dx \\ + \int_{X_{0,P}(\tilde{t})}^{1} |\rho_0^n(\tilde{\beta}_P) - \rho_0^n(\beta_P)| dx + \int_{X_{0,P}(\tilde{t})}^{1} |\rho_0^n(\beta_P) - \rho_0(\beta_P)| dx \\ \le 2 \|\rho_0^n - \rho_0\|_{L^1(0,1)} + \|(\rho_0^n)'\|_{\infty, [0,1]} \int_{X_{0,P}(\tilde{t})}^{1} |\tilde{\beta}_P - \beta_P| dx. $$

By the continuity of $\beta_P$ associated with the dominated convergence Theorem,

$$ \int_{X_{0,P}(\tilde{t})}^{1} |\tilde{\beta}_P - \beta_P| dx \xrightarrow{\tilde{t} \to t^+} 0. \text{ So by taking } n \text{ large enough and then } \tilde{t} - t \text{ small} $$
enough, the right-hand side of the last inequality can be smaller than any given positive constant, so:

$$ \lim_{\tilde{t} \to t^+} \|\rho(\tilde{t}, \cdot) - \rho(t, \cdot)\|_{L^1(X_{0,P}(\tilde{t}),1)} = 0. $$

Hence, as $\rho \in L^\infty((0, \delta) \times (0, 1))$, we can write:

$$ \int_{X_{0,P}(t)}^{X_{0,P}(\tilde{t})} |\rho(\tilde{t}, x) - \rho(t, x)| \, dx \le 2 \|\rho\|_{\infty, (0,\delta) \times (0,1)} |X_{0,P}(\tilde{t}) - X_{0,P}(t)|. $$

Now,

$$ |X_{0,P}(\tilde{t}) - X_{0,P}(t)| = \left| \int_t^{\tilde{t}} \Lambda(P(\theta), X_{0,P}(\theta)) d\theta \right| \le |\tilde{t} - t| \|\Lambda\|_{\infty, [0, M]^N \times [0,1]}. $$

Therefore,

$$ \lim_{\tilde{t} \to t^+} \|\rho(\tilde{t}, \cdot) - \rho(t, \cdot)\|_{L^1(X_{0,P}(t), X_{0,P}(\tilde{t}))} = 0. $$

Hence, by linearity of the integral:

$$ \lim_{\tilde{t} \to t^+} \|\rho(\tilde{t}, \cdot) - \rho(t, \cdot)\|_{L^1(0,1)} = 0. $$

To end this proof, let us show that $\rho$ verifies the weak formulation Eq. (2.1).

By definition of the mapping $F$ (see Eq. (2.7)) and the candidate $\rho$, we have directly that for all $i \in \{1, \dots, N\}$ and $t \in [0, \delta]$, $P_i(t) = \int_0^1 \rho(t, x)\omega_i(x)dx$.

Let $\tau \in [0, T]$ and $\varphi \in C^1([0, \tau] \times [0, 1])$ such that $\forall x \in [0, 1]$, $\varphi(\tau, x) = 0$, and $\forall t \in [0, \tau]$, $\varphi(t, 1) = 0$. We define:

$$ \begin{aligned} \star &:= \int_0^\tau \int_0^1 \rho(t, x) [\partial_t \varphi(t, x) + \Lambda(P(t), x) \partial_x \varphi(t, x) - \mu(P(t), x) \varphi(t, x)] \, dxdt \\ &= \int_0^\tau \int_0^{X_{0,P}(t)} \rho(t, x) [\partial_t \varphi(t, x) + \Lambda(P(t), x) \partial_x \varphi(t, x) - \mu(P(t), x) \varphi(t, x)] \, dxdt \\ &+ \int_0^\tau \int_{X_{0,P}(t)}^1 \rho(t, y) [\partial_t \varphi(t, y) + \Lambda(P(t), y) \partial_x \varphi(t, y) - \mu(P(t), y) \varphi(t, y)] \, dydt. \end{aligned} $$

By using the explicit formula of $\rho$, and then doing the change of variables $u(x) = \alpha_P(t, x)$ and $x(y) = \beta_P(t, y)$, we obtain:

$$ \begin{aligned} \star &= \int_0^\tau \int_0^t r(u, P(u)) \mathcal{F}_P(t; u, 0) [\partial_t \varphi + \Lambda \partial_x \varphi - \mu \varphi] (t, X_P(t; u, 0)) du dt \\ &+ \int_0^\tau \int_0^{\beta_P(t, 1)} \rho_0(x) \mathcal{F}_P(t; 0, x) [\partial_t \varphi + \Lambda \partial_x \varphi - \mu \varphi] (t, X_P(t; 0, x)) dx dt \\ &= \int_0^\tau \int_0^t r(u, P(u)) \frac{d[\mathcal{F}_P(t; u, 0) \varphi(t, X_P(t; u, 0)]}{dt} du dt \\ &+ \int_0^\tau \int_0^{\beta_P(t, 1)} \rho_0(x) \frac{d[\mathcal{F}_P(t; 0, x) \varphi(t, X_P(t; 0, x)]}{dt} dx dt. \end{aligned} $$
And by Fubini Theorem:

$$
\begin{aligned}
\star = & \int_0^\tau r(u, P(u)) \int_u^\tau \frac{d[\mathcal{F}_P(t; u, 0)\varphi(t, X_P(t; u, 0))]}{dt} dt du \\
& + \int_0^{\beta_P(\tau, 1)} \rho_0(x) \int_0^\tau \frac{d[\mathcal{F}_P(t; 0, x)\varphi(t, X_P(t; 0, x))]}{dt} dt dx \\
& + \int_{\beta_P(\tau, 1)}^1 \rho_0(x) \int_0^{(\beta_P(\cdot, 1))^{-1}(x)} \frac{d[\mathcal{F}_P(t; 0, x)\varphi(t, X_P(t; 0, x))]}{dt} dt dx.
\end{aligned}
$$

Finally, as $\forall x \in [0, 1]$, $\varphi(\tau, x) = 0$, and $\forall t \in [0, \tau]$, $\varphi(t, 1) = 0$, we get:

$$
\begin{aligned}
\star = & - \int_{0}^{\tau} r(u, P(u))\varphi(u, 0)du - \int_{0}^{\beta_P(\tau, 1)} \rho_0(x)\varphi(0, x)dx - \int_{\beta_P(\tau, 1)}^{1} \rho_0(x)\varphi(0, x)dx \\
& = - \int_{0}^{\tau} r(u, P(u))\varphi(u, 0)du - \int_{0}^{1} \rho_0(x)\varphi(0, x)dx,
\end{aligned}
$$

so $\rho$ verifies the weak formulation Eq. (2.1).

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

End of the proof of Lemma 2.4. First, let us check that, for all $P \in \Omega_{\delta,M}$, $F(P)$ is included in $C^0([0, \delta])^N$. Let $i \in \{1, \dots, N\}$, then:

$$
F_i(P)(t) = \int_0^t r(u, P(u)) \mathcal{F}_P(t; u, 0) \omega_i(X_P(t; u, 0)) du \\
\qquad + \int_0^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta} \rho_0(x) \mathcal{F}_P(t; 0, x) \omega_i(X_P(t; 0, x)) dx.
$$

For sake of readability, we note $\phi(t, u) = r(u, P(u))\mathcal{F}_P(t; u, 0)\omega_i(X_P(t; u, 0))$ and $\psi(x, u) = \rho_0(x)\mathcal{F}_P(t; 0, x)\omega_i(X_P(t; 0, x))$. From regularity assumptions (H2), (H3), (H4), (H6), we can see that $t \to \phi(t, u)$ and $t \to \psi(t, x)$ are continuous for both ev-
ery $u \in [0, \delta]$ and $x \in [0, 1 - \int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta]$. To check the continuity of
$F(P)$, we will use the sequential characterization of continuity. We have to consider
the left continuity separately from the right continuity, in order to respect the domain
of definition of the characteristic curves. First, let $t_n \xrightarrow{n \to \infty} t^-$, then:

$$
\begin{aligned}
\left| \int_0^{t_n} \phi(t_n, u) du - \int_0^t \phi(t, u) du \right| &= \left| \int_0^{t_n} \phi(t_n, u) - \phi(t, u) du + \int_0^{t_n} \phi(t, u) du \right| \\
&\le \int_0^t 1_{[0, t_n]}(u) |\phi(t_n, u) - \phi(t, u)| du \\
&\quad + |t_n - t| \|r\|_{\infty, [0, T] \times [0, M]^N} \|\omega\|_{\infty},
\end{aligned}
$$

and $\phi(t_n, u) \le \|r\|_{\infty, [0,T] \times [0,M]^N} \|\omega\|_{\infty}$, such that we can apply the dominated convergence theorem to the first right-hand side term, which leads to:

$$
\lim_{t_n \to t^-} \int_0^{t_n} \phi(t_n, u) du = \int_0^t \phi(t, u) du.
$$
Furthermore,

$$
\begin{aligned}
& \left| \int_0^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta} \psi(t_n, u)dx - \int_0^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta} \psi(t, u)dx \right| \\
& \le \int_0^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta} |\psi(t_n, u) - \psi(t, u)|dx + \int_{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta}^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta)} \psi(t_n, u)dx \\
& \le \int_0^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta} |\psi(t_n, x) - \psi(t, x)|dx + |t_n - t| \|\rho_0\|_{\infty} \|\omega\|_{\infty} \|\Lambda\|_{\infty, [0, M]^N \times [0, 1]},
\end{aligned}
$$

and, as $\psi(t_n, x) \le \|\rho_0\|_{\infty} \|\omega\|_{\infty}$, we can again apply the dominated convergence theorem to the first right-hand side term, which leads to:

$$ \lim_{n \to +\infty} \int_0^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta} \psi(t_n, x)dx = \int_0^{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta} \psi(t, u)dx. $$

Therefore $\lim_{n \to +\infty} F_i(P)(t_n) = F_i(P)(t)$. For $t_n \to t^+$, the arguments are completely similar.

Now, let us take a closer look at the calculations to show that $F$ is a contraction map, starting from inequality Eq. (2.8). For the first r.h.s. term, using in particular (H4) and (H6), we have:

$$
\begin{aligned}
& \int_0^t |r(u, P(u)) - r(u, Q(u))| \mathcal{F}_P(t; u, 0) \omega_i(X_P(t; u, 0)) du \\
& \le \|\omega\|_{\infty} \int_0^t L_{r, [0, T] \times [0, M]^N} \|P(u) - Q(u)\|_{\mathbb{R}^N} du \\
& \le \delta \|\omega\|_{\infty} L_{r, [0, T] \times [0, M]^N} \|P - Q\|_{\Omega_{\delta, M}}.
\end{aligned}
$$

For the second r.h.s. term, using (H3), (H4), (H6):

$$
\begin{aligned}
& \int_0^t r(u, Q(u)) |\mathcal{F}_P(t; u, 0) - \mathcal{F}_Q(t; u, 0)| \omega_i(X_P(t; u, 0)) du \\
& \le \|\omega\|_{\infty} \|r\|_{\infty, [0, T] \times [0, M]^N} \int_0^t |\mathcal{F}_P(t; u, 0) - \mathcal{F}_Q(t; u, 0)| du.
\end{aligned}
$$

And for all $0 \le u \le t \le \delta$,

$$
\begin{aligned}
|\mathcal{F}_P(t; u, 0) - \mathcal{F}_Q(t; u, 0)| &= \left| e^{-\int_u^t \mu(P(\theta), X_P(\theta; u, 0))d\theta} - e^{-\int_u^t \mu(Q(\theta), X_P(\theta; u, 0))d\theta} \right| \\
&\le \left| \int_u^t \mu(P(\theta), X_P(\theta; u, 0))d\theta - \int_u^t \mu(Q(\theta), X_Q(\theta; u, 0))d\theta \right| \\
&\le \int_u^t L_{\mu, M} (\|P(\theta) - Q(\theta)\|_{\mathbb{R}^N} + |X_P(\theta; u, 0) - X_Q(\theta; u, 0)|) d\theta \\
&\le \delta L_{\mu, M} \|P - Q\|_{\Omega_{\delta, M}} + \delta \|X_P(\cdot, 0; x) - X_P(\cdot, 0; x)\|_{\infty, [0, \delta]},
\end{aligned}
$$

with $L_{\mu, M} := L_{\mu, [0, M]^N \times [0, 1]}$. By Lemma 2.3,

$$ |\mathcal{F}_P(t; u, 0) - \mathcal{F}_Q(t; u, 0)| \le \delta L_{\mu, M} \|P - Q\|_{\Omega_{\delta, M}} + 2\delta^2 L_{\Lambda, M} \|P - Q\|_{\Omega_{\delta, M}}, $$
so,

$$ \int_0^t r(u, Q(u)) |\mathcal{F}_P(t; u, 0) - \mathcal{F}_Q(t; u, 0)| \omega_i(X_P(t; u, 0)) du \\ \le \|\omega\|_{\infty} \|r\|_{\infty, [0,T] \times [0,M]^N} (\delta L_{\mu,M} + 2\delta^2 L_{\Lambda,M}) \|P - Q\|_{\Omega_{\delta,M}}. $$

For the third r.h.s. term, with (H4), (H6) and Lemma 2.3:

$$ \int_0^t r(u, Q(u)) \mathcal{F}_Q(t; u, 0) |\omega_i(X_P(t; u, 0)) - \omega_i(X_Q(t; u, 0))| du \\ \le \|r\|_{\infty, [0,T] \times [0,M]^N} \int_0^t |\omega_i(X_P(t; u, 0)) - \omega_i(X_Q(t; u, 0))| du \\ \le \|r\|_{\infty, [0,T] \times [0,M]^N} L_{\omega} 2\delta^2 L_{\Lambda,M} \|P - Q\|_{\Omega_{\delta,M}}. $$

We have already shown that

$$ \left| \int_{1-\int_0^t \Lambda((P(\theta), X_P(\theta; t, 1))d\theta}^{1-\int_0^t \Lambda((Q(\theta), X_Q(\theta; t, 1))d\theta)} \rho_0(x) \mathcal{F}_P(t; 0, x) \omega_i(X_P(t; 0, x)) dx \right| \\ \le \delta \|\rho_0\|_{\infty} \|\omega\|_{\infty} L_{\Lambda,M} (1 + 2\delta L_{\Lambda,M}) \|P - Q\|_{\Omega_{\delta,M}}. $$

Furthermore, in a similar way to what we did for the second r.h.s. term:

$$ \int_0^{1-\int_0^t \Lambda((Q(\theta), X_Q(\theta; t, 1))d\theta)} \rho_0(x) |\mathcal{F}_P(t; 0, x) - \mathcal{F}_Q(t; 0, x)| \omega_i(X_P(t; 0, x)) dx \\ \le \|\rho_0\|_{\infty} \|\omega\|_{\infty} (\delta L_{\mu,M} + 2\delta^2 L_{\mu,M} L_{\Lambda,M}) \|P - Q\|_{\Omega_{\delta,M}}. $$

Finally, as for the third r.h.s. term:

$$ \int_0^{1-\int_0^t \Lambda((Q(\theta), X_Q(\theta; t, 1))d\theta)} \rho_0(x) \mathcal{F}_Q(t; 0, x) |\omega_i(X_P(t; 0, x)) - \omega_i(X_Q(t; 0, x))| dx \\ \le 2\delta \|\rho_0\|_{\infty} L_{\omega} L_{\Lambda,M} \|P - Q\|_{\Omega_{\delta,M}}. $$

Hence, by adding up all these inequalities:

$$ \begin{aligned} \|F(P) - F(Q)\|_{\Omega_{\delta,M}} &\le \delta N \left[ \|\omega\|_{\infty} L_{r, [0,T] \times [0,M]^N} \right. \\ &\quad + \|\omega\|_{\infty} \|r\|_{\infty, [0,T] \times [0,M]^N} (L_{\mu,M} + 2TL_{\Lambda,M}) \\ &\quad + \|r\|_{\infty, [0,T] \times [0,M]^N} L_{\omega} 2TL_{\Lambda,M} \\ &\quad + \|\rho_0\|_{\infty} \|\omega\|_{\infty} L_{\Lambda,M} (1 + 2TL_{\Lambda,M}) \\ &\quad + \|\rho_0\|_{\infty} \|\omega\|_{\infty} (L_{\mu,M} + 2TL_{\mu,M}L_{\Lambda,M}) \\ &\quad + 2 \|\rho_0\|_{\infty} L_{\omega} L_{\Lambda,M} \|P - Q\|_{\Omega_{\delta,M}}. \end{aligned} \quad (A.5) $$

Therefore, by taking $\delta$ small enough, $F$ is a contraction map. □

**Appendix B. Sketch of the proof of existence and uniqueness for $L^1$ initial condition.** We now aim to briefly explain how, under additional assumptions (H8)-(H10), it can be shown that a unique weak solution to Eq. (3.1) in $C^0([0, T]; L^1(0, 1))$ exists for any $L^1$ initial condition (instead of $L^\infty$ in Theorem 2.2).
First, let us note that, thanks to Proposition 3.2, this amounts to proving the existence of a unique solution to (3.4).

By analogy with subsection 2.2, we define, for any $(s, t, x) \in [0, T] \times [0, T] \times (0, 1)$,

$$X(s; t, x) = s - t + x,$$

and

$$\mathcal{F}(s; t, x) = e^{-\int_{t}^{s} m(X(\theta; t, x)) d\theta}.$$

As for the case with $L^\infty$ initial condition, the key point is to show the following Lemma B.1, which is analogous to Lemma 2.4.

**LEMMA B.1.** For some $\delta > 0$ and $M > 0$, the map $F$ defined by:

$$ (B.1) \qquad F: \Omega_{\delta,M} \to \Omega_{\delta,M} \\ P \mapsto (F_1(P), \dots, F_N(P)) $$

with, for all $i \in \{1, \dots, N\}$, $F_i$ defined by

$$ (B.2) \qquad F_i(P)(t) := \int_0^t b(P(u)) \mathcal{F}(t; u, 0) \alpha_i(X(t; u, 0)) du \\ \qquad \qquad + \int_0^{\tau-t} v_0(x) \mathcal{F}(t; 0, x) \alpha_i(X(t; 0, x)) dx, $$

is a contraction mapping.

**Proof.** First, by direct upper bounding of Eq. (B.2), we can see (with boundedness assumption (H10)) that, for all $i \in \{1, \dots, N\}$, and $0 < \delta < T$:

$$\|F_i(P)\|_{C^0([0,\delta])} \leq T e^{-Tm} \|\alpha_i\|_{\infty} \|b\|_{\infty, [0, +\infty)^N} + e^{-Tm} \|\alpha_i\|_{\infty} \|v_0\|_{L^1(0,\tau)}$$

Then, if we define $M = T e^{-Tm} \|\alpha_i\|_{\infty} \|b\|_{\infty, [0, +\infty)^N} + e^{-Tm} \|\alpha_i\|_{\infty} \|v_0\|_{L^1(0,\tau)}$, we have $F(\Omega_{\delta,M}) \subset \Omega_{\delta,M}$.

Now, let us show that $F$ is a contraction mapping. First, we can easily see from (H10) that $b$ is globally Lipschitz continuous (with Lipschitz constant denoted by $L_b$). Therefore, let $P, Q \in \Omega_{\delta,M}$, $i \in \{1, \dots, N\}$, and $t \in [0, \delta]$:

$$ \begin{aligned} |F_i(P)(t) - F_i(Q)(t)| &\leq \int_0^t |b(P(u)) - b(Q(u))| \mathcal{F}(t; u, 0) \alpha_i(X(t; u, 0)) du \\ &\leq L_b \int_0^t \|P(u) - Q(u)\| \mathcal{F}(t; u, 0) \alpha_i(X(t; u, 0)) du \\ &\leq L_b \|\alpha_i\|_{\infty} e^{-Tm} \int_0^t \|P(u) - Q(u)\| du \\ &\leq \delta L_b \|\alpha_i\|_{\infty} e^{-Tm} \|P(u) - Q(u)\| . \end{aligned} $$

Hence, if we take $\delta$ small enough, $F$ is a contraction map. $\square$

Let us notice that $\delta$ does not depend on $\|\rho_0\|_\infty$, so we can extend the local solution to a global one by iterating without needing an upper bounding on $\|\rho\|_\infty$.
## Appendix C. Time and Space Scaling Transformation.

DEFINITION C.1 ($L^1$ - stability). A stationary solution $\bar{\rho}$ of Eq. (3.1) is locally stable if

$$ \forall \epsilon > 0, \exists \delta_{\epsilon} > 0 : \forall \rho_{0} \in L^{1}(0,1), (\|\bar{\rho} - \rho_{0}\|_{L^{1}} \le \delta_{\epsilon}) \Rightarrow (\forall t \ge 0, \|\rho(t, \cdot) - \bar{\rho}\|_{L^{1}} \le \epsilon). $$

Furthermore, $\bar{\rho}$ is locally asymptotically stable if

$$ \exists \delta > 0 : \forall \rho_0 \in L^1(0,1), (\|\bar{\rho} - \rho_0\|_{L^1} \le \delta) \Rightarrow \left( \|\rho(t, \cdot) - \bar{\rho}\|_{L^1} \xrightarrow{t \to +\infty} 0 \right). $$

If $\bar{\rho}$ is not locally stable, it is said unstable.

### Proof of Proposition 3.2.
Let $\rho$ be a weak solution to Eq. (3.1) up to time $T$. By regularity assumption (H2) and positivity assumption (H8), $\Gamma : x \mapsto \int_0^x \frac{1}{f(y)} dy$ and $s_\rho$ defined by Eq. (3.5) are $C^1$ bijective functions respectively from $[0, 1]$ to $[0, \tau = \Gamma(1)]$ and from $[0, T]$ to $[0, s_\rho(T) := S]$. Therefore, we can define $v \in C^0([0, S], L^1(0, \tau))$ such that $v(s, y) = \rho(s_\rho^{-1}(s), \Gamma^{-1}(y))$. Furthermore, it is easy to check that $s_\rho^{-1} = t_v$, with $t_v$ defined by Eq. (3.6).

Now, let $\sigma \in [0, S]$, and $\phi \in C^1([0, \sigma] \times [0, \tau])$ such that $\phi(\sigma, y) = 0$ for all $y \in [0, \tau]$ and $\phi(s, \tau) = 0$ for all $s \in [0, \sigma]$. Let us define.

$$ \begin{aligned} \diamond := \int_0^\sigma \int_0^\tau v(s, y) [\partial_s \phi(s, y) + \partial_y \phi(s, y) - m(y)\phi(s, y)] dy ds \\ + \int_0^{t_v(\sigma)} \frac{1}{f(0)} b(V(s)) \phi(s, 0) ds + \int_0^\tau \rho_0(\Gamma^{-1}(y)) \phi(0, y) dy. \end{aligned} $$

With the changes of variables $t = t_v(s)$, and $x = \Gamma^{-1}(y)$, we define $\psi(t, x) := \phi(s_\rho(t), \Gamma(x))$. By composition, $\psi \in C^1([0, t_v(\sigma)] \times [0, 1])$, $\partial_t \psi(t, \Gamma(x)) = g(P(t)) \partial_s \phi(s_\rho(t), \Gamma(x))$, and $\partial_x \psi(t, x) = \frac{1}{f(x)} \partial_y \phi(s_\rho(t), \Gamma(x))$. Thus,

$$ \begin{aligned} \diamond = \int_0^{t_v(\sigma)} \int_0^1 \rho(t, x) \left[ \frac{\partial_t \psi(t, x)}{g(P(t))} + f(x) \partial_x \psi(t, x) - f'(x) \psi(t, x) \right] \frac{g(P(t))}{f(x)} dx dt \\ + \int_0^\sigma \frac{1}{f(0)} r(P(t)) \psi(t, 0) dt + \int_0^1 \rho_0(x) \frac{\psi(0, x)}{f(x)} dx \\ = \int_0^{t_v(\sigma)} \int_0^1 \rho(t, x) \left[ \partial_t \left( \frac{\psi(t, x)}{f(x)} \right) + f(x) g(P(t)) \left( \frac{\partial_x \psi(t, x)}{f(x)} - \frac{f'(x)}{f^2(x)} \psi(t, x) \right) \right] dx dt \\ + \int_0^\sigma \frac{\psi(t, 0)}{f(0)} r(P(t)) dt + \int_0^1 \rho_0(x) \frac{\psi(0, x)}{f(x)} dx \\ = \int_0^{t_v(\sigma)} \int_0^1 \rho(t, x) \left[ \partial_t \left( \frac{\psi(t, x)}{f(x)} \right) + f(x) g(P(t)) \partial_x \left( \frac{\psi(t, x)}{f(x)} \right) \right] dx dt \\ + \int_0^\sigma \frac{\psi(t, 0)}{f(0)} r(P(t)) dt + \int_0^1 \rho_0(x) \frac{\psi(0, x)}{f(x)} dx. \end{aligned} $$

As $\rho$ is a weak solution to Eq. (3.1) and $(t, x) \mapsto \psi(t, x)/f(x)$ is actually a test function, $\diamond = 0$ hence $v$ is a weak solution to Eq. (3.4). The reciprocal works in the same way. □
## Appendix D. Spectral Properties of $A + DG(\bar{v})$.

*Proof details of Proposition 3.5, part 1.* Firstly, we can observe that $A$ can be decomposed as the sum of two operators $A = A_1 + A_2$, where

$$A_1 \begin{pmatrix} 0 \\ \phi \end{pmatrix} = \begin{pmatrix} -\phi(0) \\ -\partial_y \phi \end{pmatrix}, \text{ is closed (direct calculations), and}$$

$$A_2 \begin{pmatrix} 0 \\ \phi \end{pmatrix} = \begin{pmatrix} 0 \\ -m\phi \end{pmatrix} \text{ is bounded as } m \in C^0(0, \tau).$$

Then $A$ is closed as the sum of a linear closed operator and a bounded operator. Consider $\lambda > -\underline{m}$, and $(a, \psi) \in \mathbb{R} \times L^1(0, \tau)$:

$$ (\lambda \mathbb{I} - A)^{-1} \begin{pmatrix} a \\ \psi \end{pmatrix} = \begin{pmatrix} 0 \\ \phi \end{pmatrix} \Leftrightarrow \left\{ \begin{array}{l} \phi(0) = a \\ \phi' + (m + \lambda)\phi = \psi \end{array} \right. \\ \Leftrightarrow \phi(y) = ae^{-\int_0^y \lambda+m(x)dx} + \int_0^y \psi(s)e^{-\int_s^y \lambda+m(x)dx} ds, \quad y \in [0, \tau]. $$

We can compute

$$ \begin{aligned} \left\| (\lambda \mathbb{I} - A)^{-1} \begin{pmatrix} a \\ \psi \end{pmatrix} \right\|_{X} &= \left\| \begin{pmatrix} 0 \\ \phi \end{pmatrix} \right\|_{X} = |0| + \|\phi\|_{L^1(0, \tau)} \\ &\le \int_0^\tau |a| e^{-\int_0^y \lambda+m(x)dx} dy + \int_0^\tau \int_0^y |\psi(s)| e^{-\int_s^y \lambda+m(x)dx} ds dy \\ &\le |a| \int_0^\tau e^{-(\lambda+\underline{m})y} dy + \int_0^\tau |\psi(s)| \int_s^\tau e^{-\int_s^y \lambda+m(x)dx} dy ds \\ &\le |a| \frac{1 - e^{-(\lambda+\underline{m})\tau}}{\lambda + \underline{m}} + \int_0^\tau |\psi(s)| \int_s^\tau e^{-(y-s)(\lambda+\underline{m})} dy ds \\ &\le |a| \frac{1 - e^{-(\lambda+\underline{m})\tau}}{\lambda + \underline{m}} + \int_0^\tau |\psi(s)| \frac{1 - e^{-(\tau-s)(\lambda+\underline{m})}}{\lambda + \underline{m}} ds \\ &\le \frac{1}{\lambda + \underline{m}} \left( |a| + \int_0^\tau |\psi(s)| ds \right) = \frac{1}{\lambda + \underline{m}} \left\| \begin{pmatrix} a \\ \psi \end{pmatrix} \right\|_{X}. \end{aligned} $$

Then

$$ \left\| (\lambda \mathbb{I} - A)^{-1} \right\|_{\mathcal{L}(X)} \le \frac{1}{\lambda + \underline{m}}, $$

hence, for all $n \ge 1$,

$$ \left\| (\lambda \mathbb{I} - A)^{-n} \right\|_{\mathcal{L}(X)} \le \left\| (\lambda \mathbb{I} - A)^{-1} \right\|_{\mathcal{L}(X)}^n \le \frac{1}{(\lambda + \underline{m})^n}, $$

and $\lambda \in \rho(A)$.

□

The following Lemma will be useful for studying the spectrum of $A + DG(\bar{v})$.

**LEMMA D.1.** Let $\lambda \in \mathbb{C}$, $\lambda I - A - DG(\bar{v}) : D(A) \to X$ is invertible if and only if $I - DG(\bar{v}) (\lambda I - A)^{-1} : X \to X$ is invertible, and in this case,

$$ (\lambda I - A - DG(\bar{v}))^{-1} = (\lambda I - A)^{-1} (I - DG(\bar{v}) (\lambda I - A)^{-1})^{-1}. $$
The proof is straightforward. We now turn to the second part of the proof of Proposition 3.5.

*Proof details of Proposition 3.5, part 2.* It remains only to determine the spectrum of $A + DG(\bar{v})$ and its properties. Firstly, we can notice that $\sigma(A + DG(\bar{v})) = \sigma((A + DG(\bar{v}))_{\overline{D(A)}})$ (see Lemma 2.2.10. in [37]). Let $\begin{pmatrix} a \\ \phi \end{pmatrix} \in X$, then

$$ (\lambda I - A)^{-1} \begin{pmatrix} a \\ \phi \end{pmatrix} = \begin{pmatrix} 0 & 0 \\ y \mapsto ae^{-\int_0^y \lambda+m(x)dx} & + \int_0^y \phi(s)e^{-\int_s^y \lambda+m(x)dx} ds \end{pmatrix}, $$

hence

$$ \begin{aligned} & (I - DG(\bar{v}) (\lambda I - A)^{-1}) \begin{pmatrix} a \\ \phi \end{pmatrix} \\ &= \left( a - \frac{1}{f(0)} \langle \nabla_x b(\gamma, \bar{v}_\gamma) \int_0^\tau \alpha(y) \left[ ae^{-\int_0^y \lambda+m(x)dx} + \int_0^y \phi(s)e^{-\int_s^y \lambda+m(x)dx} ds \right] dy \right) \begin{pmatrix} a \\ \phi \end{pmatrix} \\ &= \left( a\Delta(\gamma, \lambda) - \frac{1}{f(0)} \langle \nabla_x b(\gamma, \bar{v}_\gamma) \int_0^\tau \alpha(y) \int_0^y \phi(s)e^{-\int_s^y \lambda+m(x)dx} ds dy \right) \end{aligned} $$

with

$$ \begin{aligned} \Delta(\gamma, \lambda) &:= 1 - \langle \nabla_x b(\gamma, \bar{v}_\gamma), \int_0^1 \frac{\omega(x)}{f(x)} e^{-\lambda \Gamma(x)} dx \rangle \\ &= 1 - \frac{1}{f(0)} \langle \nabla_x b(\gamma, \bar{v}_\gamma), \int_0^\tau \alpha(y) e^{-\int_0^y \lambda+m(x)dx} dy \rangle. \end{aligned} $$

So it is easy to see that if $\Delta(\gamma, \lambda) \neq 0$, then $I - DG(\bar{v}) (\lambda I - A)^{-1}$ is invertible, hence by Lemma D.1, $\lambda I - A - DG(\bar{v})$ is invertible, i.e. $\lambda \in \rho(A + DG(\bar{v}))$.

By contraposition,

$$ \sigma(A + DG(\bar{v})) \subset \{\lambda \in \mathbb{C} : \Delta(\gamma, \lambda) = 0\}. $$

Now, let us consider $\lambda \in \mathbb{C}$ such that $\Delta(\gamma, \lambda) = 0$. We set $\begin{pmatrix} 0 \\ \phi \end{pmatrix} \in D(A) \setminus \{0\}$ such that:

$$ \begin{pmatrix} 0 \\ \phi \end{pmatrix} := (\lambda I - A)^{-1} \begin{pmatrix} a \\ 0 \end{pmatrix} = \begin{pmatrix} 0 & 0 \\ y \mapsto ae^{-\int_0^y \lambda+m(x)dx} & \end{pmatrix} $$

with $a \neq 0$. Then

$$ \begin{aligned} (\lambda I - A - DG(\bar{v})) \begin{pmatrix} 0 \\ \phi \end{pmatrix} &= (\lambda I - A - DG(\bar{v})) (\lambda I - A)^{-1} \begin{pmatrix} a \\ 0 \end{pmatrix} \\ &= (I - DG(\bar{v}) (\lambda I - A)^{-1}) \begin{pmatrix} a \\ 0 \end{pmatrix} \\ &= \begin{pmatrix} a \Delta(\gamma, \lambda) \\ 0 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \end{aligned} $$
so $\begin{pmatrix} 0 \\ \phi \end{pmatrix}$ is an eigenvector of $A + DG(\bar{v})$ associated with the eigenvalue $\lambda$, and then

$$ \{\lambda \in \mathbb{C} : \Delta(\gamma, \lambda) = 0\} \subset \sigma_p(A + DG(\bar{v})) \subset \sigma(A + DG(\bar{v})). $$

Finally,

$$ \sigma(A + DG(\bar{v})) = \sigma_p(A + DG(\bar{v})) = \{\lambda \in \mathbb{C} : \Delta(\gamma, \lambda) = 0\}. \quad \square $$

*Proof details of Corollary 3.7.* The $C_0$-semigroup generated by $A+DG(\bar{v})$ is positive if and only if the resolvent of $R(\lambda, A+DG(\bar{v}))$ is positive for all sufficiently large $\lambda$ (see [18] Theorem VI.1.8.).

First at all, notice that

$$ \lim_{\lambda \to +\infty} \Delta(\gamma, \lambda) = 1, $$

so for any sufficiently large $\lambda$, $\Delta(\gamma, \lambda) \neq 0$ and $R(\lambda, A+DG(\bar{v}))$ is well defined.

Let us consider $\lambda \in \mathbb{C}$ such that $\Delta(\gamma, \lambda) \neq 0$. Then, thanks to Lemma D.1, we have

$$ R(\lambda, A + DG(\bar{v})) = (\lambda I - A)^{-1} (I - DG(\bar{v}) (\lambda I - A)^{-1})^{-1}. $$

Yet,

$$ (I - DG(\bar{v}) (\lambda I - A)^{-1})^{-1} \begin{pmatrix} a \\ \phi \end{pmatrix} = \begin{pmatrix} \frac{a + \eta(\phi, \lambda)}{\Delta(\gamma, \lambda)} \\ \phi \end{pmatrix} $$

where

$$ \eta(\phi, \lambda) := \frac{1}{f(0)} \langle \nabla_x, b(\gamma, \bar{v}_\gamma) \int_0^\tau \alpha(y) \int_0^y \phi(s) e^{-\int_s^y \lambda + m(x)dx} ds dy \rangle. $$

Let us suppose that $a \ge 0$ and $\phi \ge 0$. Then, it is easy to see that if Eq. (3.15) is satisfied, then $\eta(\phi, \lambda) \ge 0$ and $\frac{a+\eta(\phi,\lambda)}{\Delta(\gamma,\lambda)} \ge 0$. So $I - DG(\bar{v}) (\lambda I - A)^{-1}$ is positive for all sufficiently large $\lambda$. Furthermore, always with $a \ge 0$ and $\phi \ge 0$, we have directly

$$ (\lambda I - A)^{-1} \begin{pmatrix} a \\ \phi \end{pmatrix} = \begin{pmatrix} 0 \\ \int_0^y \lambda + m(x) dx + \int_0^y \phi(s) e^{-\int_s^y \lambda + m(x) dx} ds \end{pmatrix} \ge 0. $$

Thus $R(\lambda, A)$ is positive too. Finally, by composition of positive operators, $R(\lambda, A + DG(\bar{v}))$ is positive for any sufficiently large $\lambda$, and then the $C_0$-semigroup generated by $A + DG(\bar{v})$ is positive. $\square$

## Appendix E. Characteristic Equation of the Spawning Cycle Model.

*Proof of Proposition 4.4.* From Eq. (3.13), and Proposition 3.3, we recall that the characteristic function is given by

$$ \Delta(\lambda) = 1 - b'(\bar{V}) \int_0^1 \frac{\omega(x)}{f(x)} e^{-\lambda \Gamma(x)} dx = 1 - b'(\bar{P}) \int_0^1 \frac{\omega(x)}{f(x)} e^{-\lambda \Gamma(x)} dx. $$

If $f(x) \equiv 1$, then for all $x \in [0, 1]$, $\Gamma(x) = \int_0^x \frac{1}{f(x)} dx = x$, and $\omega(x) = 2Mx$ for all $x \in [0, 1/2]$, $\omega(x) = 2M(1-x)$ for all $x \in [1/2, 1]$. Then we have:

$$ \begin{aligned} \Delta(\lambda) &= 1 - 2Mb'(\bar{P}) \left[ \int_0^{1/2} xe^{-\lambda x} dx + \int_{1/2}^1 (1-x)e^{-\lambda x} dx \right] \\ &= 1 - 2Mb'(\bar{P}) \left[ \int_0^{1/2} xe^{-\lambda x} dx + e^{-\lambda} \int_0^{1/2} ye^{\lambda y} dy \right]. \end{aligned} $$
If $\lambda \neq 0$, by integration by parts

$$I(\lambda) := \int_{0}^{1/2} xe^{-\lambda x} dx = \frac{1}{\lambda^2} \left(1 - \frac{\lambda}{2} e^{-\lambda/2} - e^{-\lambda/2}\right)$$

and if $\lambda = 0$, then $I(0) = \int_{0}^{1/2} x dx = \frac{1}{8}$. So

$$\begin{aligned} \Delta(\lambda) &= 1 - 2Mb'(\bar{P}) [I(\lambda) + e^{-\lambda}I(-\lambda)] \\ &= \begin{cases} 1 - 2Mb'(\bar{P}) \frac{(1-e^{-\lambda/2})^2}{\lambda^2} & \text{if } \lambda \neq 0, \\ 1 - \frac{1}{2}Mb'(\bar{P}) & \text{if } \lambda = 0. \end{cases} \end{aligned}$$

Therefore the roots $\lambda$ of $\Delta$ verify

$$\frac{e^{-\frac{\lambda}{2}} - 1}{\lambda} = \pm i \sqrt{\frac{-1}{2b'(\bar{P})M}}.$$

Let us assume that $\lambda$ is a non trivial purely imaginary root, i.e. $\lambda = ix$ with $x \in \mathbb{R}^*$, and let us denote $c := \sqrt{\frac{-1}{2b'(\bar{P})M}}$. Then $x$ verifies

$$e^{-ix/2} \pm cx = 1,$$

so

$$(E.1) \qquad \cos(x/2) = 1 \mp cx,$$

$$(E.2) \qquad \sin(x/2) = 0.$$

From Eq. (E.2), we deduce that $x = 2q\pi$, $q \in \mathbb{Z}^*$, and Eq. (E.1) becomes:

$$(-1)^q - 1 = \mp 2q\pi c.$$

If $q$ is even, then $q = 0$ which is a contradiction. If $q$ is odd, i.e. $q = 2k + 1$,

$$k\pi + \frac{\pi}{2} = \pm \frac{1}{2c}.$$

Now, by derivating Eq. (4.3)

$$b'(\bar{P}) = \frac{-c_r c_\lambda n \bar{P}^{n-1}}{\left(1 + c_\lambda \bar{P}^n\right)^2}$$

and using Eq. (4.4), we find:

$$c := \sqrt{\frac{-1}{2b'(\bar{P})M}} = \frac{1}{2\sqrt{n}} \sqrt{\frac{Mc_r}{Mc_r - 2\bar{P}}}.$$

Finally, the characteristic function admits a pair of (non trivial) conjugated purely imaginary roots $(\lambda, \bar{\lambda} = -\lambda)$ given by $\lambda = 2q\pi i = (4k\pi + 2\pi) i$, $k \in \mathbb{N}^*$, if and only if $k$ verifies

$$k\pi + \frac{\pi}{2} = \frac{1}{2c} = \sqrt{n} \sqrt{1 - \frac{2 \bar{P}(c_{r,k})}{M c_{r,k}}}.$$
**Acknowledgment.** We would like to express our gratitude to Violette Thermes for her valuable insights and discussions on the mechanisms of fish oogenesis. We also thank the referees for their thorough reviews.

REFERENCES

[1] A. ACKLEH AND K. ITO, *An implicit finite difference scheme for the nonlinear size-structured population model*, Numer. Funct. Anal. Optim., 18 (1997), pp. 865-884.

[2] B. AYMARD, F. CLÉMENT, D. MONNIAUX, AND M. POSTEL, *Cell-kinetics based calibration of a multiscale model of structured cell populations in ovarian follicles*, SIAM J. Appl. Math., 76 (2016), pp. 1471-1491.

[3] G. BALLIF, F. CLÉMENT, AND R. YVINEC, *Averaging of a stochastic slow-fast model for population dynamics: application to the development of ovarian follicles*, SIAM J. Appl. Math., 82 (2022), pp. 359-380.

[4] G. BALLIF, F. CLÉMENT, AND R. YVINEC, *Nonlinear compartmental modeling to monitor ovarian follicle population dynamics on the whole lifespan*, J. Math. Biol., 89 (2024), p. 9.

[5] C. BARRIL, À. CALSINA, O. DIEKMANN, AND J. Z. FARKAS, *On the formulation of size-structured consumer resource models (with special attention for the principle of linearized stability)*, Math. Models Methods Appl. Sci., 32 (2022), pp. 1141-1191.

[6] C. BONNET, K. CHA HOUR, F. CLÉMENT, M. POSTEL, AND R. YVINEC, *Multiscale population dynamics in reproductive biology: singular perturbation reduction in deterministic and stochastic models*, ESAIM, Proc. Surv., 67 (2020), pp. 72-99.

[7] D. BREDA, O. DIEKMANN, M. GYLLENBERG, F. SCARABEL, AND R. VERMIGLIO, *Pseudospectral Discretization of Nonlinear Delay Equations: New Prospects for Numerical Bifurcation Analysis*, SIAM J. Appl. Dyn. Syst., 15 (2016), pp. 1-23.

[8] À. CALSINA, O. DIEKMANN, AND J. Z. FARKAS, *Structured populations with distributed recruitment: from PDE to delay formulation*, Math. Methods Appl. Sci., 39 (2016), pp. 5175-5191.

[9] À. CALSINA AND J. SALDANA, *A model of physiologically structured population dynamics with a nonlinear individual growth rate*, J. Math. Biol., 33 (1995).

[10] À. CALSINA AND M. SANCHÓN, *Stability and instability of equilibria of an equation of size structured population dynamics*, J. Math. Anal. Appl., 286 (2003), pp. 435-452.

[11] J.-M. CORON, M. KAWSKI, AND Z. WANG, *Analysis of a conservation law modeling a highly re-entrant manufacturing system*, Discrete Contin. Dyn. Syst., Ser. B, 14 (2010), pp. 1337-1359.

[12] J.-M. CORON AND Z. WANG, *Controllability for a scalar conservation law with nonlocal velocity*, J. Differ. Equations, 252 (2012), pp. 181-201.

[13] J. M. CUSHING, *An introduction to structured population dynamics*, SIAM, 1998.

[14] O. DIEKMANN, P. GETTO, AND M. GYLLENBERG, *Stability and bifurcation analysis of Volterra functional equations in the light of suns and stars*, SIAM Journal on Mathematical Analysis, 39 (2008), pp. 1023-1069.

[15] O. DIEKMANN AND M. GYLLENBERG, *Equations with infinite delay: Blending the abstract and the concrete*, J. Differ. Equations, 252 (2012), pp. 819-851.

[16] M. DOUMIC, A. MARCINIAK-CZOCHRA, B. PERTHAME, AND J. P. ZUBELLI, *A structured population model of cell differentiation*, SIAM J. Appl. Math., 71 (2011), pp. 1918-1940.

[17] A. DUCROT, Z. LIU, AND P. MAGAL, *Essential growth rate for bounded linear perturbation of non-densely defined Cauchy problems*, J. Math. Anal. Appl., 341 (2008), pp. 501-518.

[18] K.-J. ENGEL AND R. NAGEL, *One-Parameter Semigroups for Linear Evolution Equations*, Springer-Verlag New York, 2000.

[19] J. Z. FARKAS, *Structured populations: The stabilizing effect of the inflow of newborns from an external source and the net growth rate*, Appl. Math. Comput., 199 (2008), pp. 547-558.

[20] J. Z. FARKAS, *Size-structured populations: immigration, (bi)stability and the net growth rate*, J. Appl. Math. Comput., 35 (2011), pp. 617-633.

[21] J. Z. FARKAS AND T. HAGEN, *Linear stability and positivity results for a generalized size-structured daphnia model with inflow*, Appl. Anal., 86 (2007), pp. 1087-1103.

[22] E. FRANCO, O. DIEKMANN, AND M. GYLLENBERG, *Modelling physiologically structured populations: renewal equations and partial differential equations*, J. Evol. Equ., 23 (2023), p. 46.

[23] P. GETTO AND M. WAURICK, *A differential equation with state-dependent delay from cell population biology*, J. Differ. Equations, 260 (2016), pp. 6176-6200.
[24] A. GRABOSCH AND H. J. A. M. HEIJMANS, *Cauchy problems with state-dependent time evolution*, Japan J. Appl. Math., 7 (1990), pp. 433-457.

[25] A. GRABOSCH AND H. J. A. M. HEIJMANS, *Production, development, and maturation of red blood cells*, in Mathematical population dynamics, O. Arino, D. E. Axelrod, and M. Kimmel, eds., CRC Press, 1991, pp. 189-210.

[26] N. D. HAYES, *Roots of the transcendental equation associated with a certain difference-differential equation*, J. Lond. Math. Soc., s1-25 (1950), pp. 226-232.

[27] M. IANNELLI, *Mathematical theory of age-structured population dynamics*, Giardini editori e stampatori, 1995.

[28] T. IWAMATSU, T. OHTA, E. OSHIMA, AND N. SAKAI, *Oogenesis in the medaka Oryzias latipes: stages of oocyte development*, Zool. Sci., 5 (1988), pp. 353-373.

[29] N. KATO, *A principle of linearized stability for nonlinear evolution equations*, Trans. Am. Math. Soc., 347 (1995), pp. 2851-2868.

[30] N. KATO, *Positive global solutions for a general model of size-dependent population dynamics*, Abstr. Appl. Anal., 5 (2000), pp. 191-206.

[31] N. KATO, *A general model of size-dependent population dynamics with nonlinear growth rate*, J. Math. Anal. Appl., 297 (2004), pp. 234-256.

[32] K. A. KIDD, P. J. BLANCHFIELD, K. H. MILLS, V. P. PALACE, R. E. EVANS, J. M. LAZORCHAK, AND R. W. FLICK, *Collapse of a fish population after exposure to a synthetic estrogen*, PNAS, 104 (2007), pp. 8897-8901.

[33] Q. LI AND P. SHANG, *Controllability for a highly re-entrant manufacturing system with local and nonlocal velocity*, Eur. J. Control, 67 (2022), p. 100716.

[34] Z. LIU, P. MAGAL, AND S. RUAN, *Hopf bifurcation for non-densely defined Cauchy problems*, Z. Angew. Math. Phys., 62 (2011), pp. 191-222.

[35] M. C. MACKEY, M. TYRAN-KAMIŃSKA, AND R. YVINEC, *Molecular distributions in gene regulatory dynamics*, J. Theor. Biol., 274 (2011), pp. 84-96.

[36] P. MAGAL AND S. RUAN, *Center manifolds for semilinear equations with non-dense domain and applications to hopf bifurcation in age structured models*, Mem. Am. Math. Soc., 202 (2009), pp. 1-80.

[37] P. MAGAL AND S. RUAN, *Theory and applications of abstract semilinear Cauchy problems*, vol. 201 of Appl. Math. Sci., Springer International Publishing, 2018.

[38] J. A. METZ AND O. DIEKMANN, *The dynamics of physiologically structured populations*, Springer-Verlag Berlin, 1986.

[39] A. RODRÍGUEZ-MARÍ, Y.-L. YAN, R. A. BREMILLER, C. WILSON, C. CAÑESTRO, AND J. H. POSTLETHWAIT, *Characterization and expression pattern of zebrafish anti-Müllerian hormone (amh) relative to sox9a, sox9b, and cyp19a1a, during gonad development*, Gene Expr. Patterns, 5 (2005), pp. 655-667.

[40] F. SCARABEL, D. BREDA, O. DIEKMANN, M. GYLLENBERG, AND R. VERMIGLIO, *Numerical Bifurcation Analysis of Physiologically Structured Population Models via Pseudospectral Approximation*, Vietnam J. Math., 49 (2021), pp. 37-67.

[41] K. SELMAN, R. A. WALLACE, A. SARKA, AND X. QI, *Stages of oocyte development in the zebrafish, Brachydanio rerio*, J. Morphol., 218 (1993), pp. 203-224.

[42] P. SHANG AND Z. WANG, *Analysis and control of a scalar conservation law modeling a highly re-entrant manufacturing system*, J. Differ. Equations, 250 (2011), pp. 949-982.

[43] H. R. THIEME, *Well-posedness of physiologically structured population models for Daphnia magna: How biological concepts can benefit by abstract mathematical analysis*, J. Math. Biol., 26 (1988), pp. 299-317.

[44] Y. TIAN, I. AL-DARABSAH, AND Y. YUAN, *Global dynamics in sea lice model with stage structure*, Nonlinear Anal., Real World Appl., 44 (2018), pp. 283-304.

[45] S. L. TUCKER AND S. O. ZIMMERMAN, *A Nonlinear Model of Population Dynamics Containing an Arbitrary Number of Continuous Structure Variables*, SIAM J. Appl. Math., 48 (1988), pp. 549-591.

[46] R. A. WALLACE AND K. SELMAN, *Cellular and dynamic aspects of oocyte growth in teleosts.*, Am. Zool., 21 (1981), pp. 325-343.

[47] G. F. WEBB, *Theory of nonlinear age-dependent population dynamics*, Marcel Dekker, 1985.

[48] Y.-L. YAN, P. BATZEL, T. TITUS, J. SYDES, T. DESVIGNES, R. BREMILLER, B. DRAPER, AND J. H. POSTLETHWAIT, *A hormone that lost its receptor: anti-müllerian hormone (AMH) in zebrafish gonad development and sex determination*, Genetics, 213 (2019), pp. 529-553.

[49] Z. ZHANG, K. WU, Z. REN, AND W. GE, *Genetic evidence for Amh modulation of gonadotropin actions to control gonadal homeostasis and gametogenesis in zebrafish and its noncanonical signalling through Bmpr2a receptor*, Development, (2020), p. dev.189811.

[50] L.-Y. ZHOU, D.-S. WANG, Y. SHIBATA, B. PAUL-PRASANTH, A. SUZUKI, AND Y. NAGAHAMA,
*Characterization, expression and transcriptional regulation of P450c17-I and -II in the
medaka, *Oryzias latipes*, Biochem. Biophys. Res. Commun., 362 (2007), pp. 619-625.*
# SUPPLEMENTARY MATERIALS: BIFURCATION ANALYSIS OF A SIZE-STRUCTURED POPULATION MODEL: APPLICATION TO OOCYTE DYNAMICS AND OVARIAN CYCLE

FRÉDÉRIQUE CLÉMENT*, LOUIS FOSTIER††, AND ROMAIN YVINEC‡‡

**SM1. Numerical Scheme.** Numerical simulations are performed with an implicit finite difference scheme inspired from [SM1], where a convergence proof for the scheme is provided. However, it is important to note that the boundary condition used in the current work is slightly different from that in [SM1].

We use a uniform discretization in both space and time, over the domain $[0, 1] \times [0, T]$, with $T = M\Delta_t$ and $1 = K\Delta_x$. The discrete approximation of the solution is denoted as $\rho_k^n \approx \rho(t_m, x_k)$ for $0 \le m \le M$ and $1 \le k \le K$. The scheme is initialized as follows:

$$ \rho_k^0 = \frac{1}{\Delta_x} \int_{x_{k-1}}^{x_k} \rho_0(x) dx, \quad k \in \{1, \dots, K\}, $$

with $\rho_0^0 = \rho_1^0$. Moreover, for all $m \in \{1, \dots, M-1\}$:

$$ \left\{ \begin{array}{l} \lambda_0^m \rho_0^{m+1} = r^m, \\ \rho_k^{m+1} = \rho_k^m - \frac{\Delta_t}{\Delta_x} \left[ \lambda_k^m \rho_k^{m+1} - \lambda_{k-1}^m \rho_{k-1}^{m+1} \right] - \Delta_t \mu_k^m \rho_k^{m+1}, \quad k \in \{1, \dots, K\}. \end{array} \right. $$

The discretized rates are defined by:

$$ \lambda_k^m = \lambda(P^m, x_k), \quad \mu_k^m = \mu(P^m, x_k), \quad r^m = r(t_m, P^m), $$

with

$$ P_j^m = \Delta_x \sum_{k=1}^{K} \omega_j(x_k) \rho_k^m, \quad \forall j \in \{1, \dots, N\}. $$

## REFERENCES

[SM1] A. ACKLEH AND K. ITO, An implicit finite difference scheme for the nonlinear size-structured population model, Numer. Funct. Anal. Optim., 18 (1997), pp. 865–884.

<footnote>* Université Paris-Saclay, Inria, Centre Inria de Saclay, 91120, Palaiseau, France (fred-erique.clement@inria.fr).</footnote>

<footnote>† PRC, INRAE, CNRS, Université de Tours, 37380 Nouzilly, France (louis.fostier@inria.fr, ro-main.yvinec@inria.fr).</footnote>