Benjamin Aymard

Frédérique Clément

Marie Postel

Adaptive mesh refinement strategy for a non conservative transport problem *

Long time simulations of transport equations raise computational challenges since they require both a large domain of calculation and sufficient accuracy. It is therefore advantageous, in terms of computational costs, to use a time varying adaptive mesh, with small cells in the region of interest and coarser cells where the solution is smooth. Biological models involving cell dynamics fall for instance within this framework and are often non conservative to account for cell division. In that case the threshold controlling the spatial adaptivity may have to be time-dependent in order to keep up with the progression of the solution. In this article we tackle the difficulties arising when applying a Multiresolution method to a transport equation with discontinuous fluxes modeling localized mitosis. The analysis of the numerical method is performed on a simplified model and numerical scheme. An original threshold strategy is proposed and validated thanks to extensive numerical tests. It is then applied to a biological model in both cases of distributed and localized mitosis.

Introduction

When studying phenomena over long periods in the context of transport equations, one has to make simulations on large computational domains so as to follow the interesting part of the solution throughout the duration of the experiment. On the other hand it can happen that the solutions of the problem occupy only a bounded and time varying portion of the domain. Such set-ups can be very costly to simulate numerically with uniform grids. Timeadaptive mesh refinement methods enable one to minimize the size of the mesh, hence both the memory requirement and CPU time, by restricting the use of fine discretization to the areas of interest, or more generally to the regions where the solution varies rapidly. Among these methods, the Multiresolution (MR) method is specially interesting for hyperbolic systems of PDEs. Firstly, the underlying wavelet theory enables one to control the numerical error induced by the adaptive scheme with a small parameter ε that is explicitly used as a threshold in the algorithm. Secondly, MR provides one with a strategy to monitor the evolution in time of the adaptive grid, relying on the hyperbolicity of the equations to predict the displacement of the singularities. We refer the reader to [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF][START_REF] Müller | Adaptive multiscale schemes for conservation laws[END_REF] for a thorough analysis of the MR method applied to a scalar hyperbolic equation and to [START_REF] Cohen | Adaptive multiresolution for finite volume solutions of gas dynamics[END_REF][START_REF] Coquel | Local time stepping applied to implicit-explicit methods for hyperbolic systems[END_REF] and references therein for applications to more realistic models and set-ups. Although originally designed for and theoretically studied in the case of hyperbolic equations, Multiresolution is currently used to simulate reaction-diffusion and parabolic equations with interesting computational performances, as illustrated in [START_REF] Coquel | Multiresolution and Adaptive Methods for Convection-Dominated Problems[END_REF][START_REF] Louvet | Summer school on multiresolution and adaptive mesh refinement methods[END_REF].

In the current work we deal with non conservative transport equations frequently arising in biological models involving cell dynamics. In these models, which are most of the time 1D, the space variable plays the role of age and marks the progression along the cell cycle. Mitosis is the endpoint of the cell cycle, that gives birth to two daughter cells from one mother cell. From the modeling viewpoint, mitosis can be distributed over the cell cycle thanks to a positive linear source term, amounting to doubling the cell density in a time equivalent to the duration of the cell cycle (see for instance [START_REF] Perthame | Transport Equations In Biology[END_REF]). Another approach treats the mitosis phenomenon more directly, as a localized event at the end of each cell cycle, corresponding to a doubling transmission condition which must be satisfied by the cell density at the interface between cell cycles (see for instance [START_REF] Shang | Cauchy problem for multiscale conservation laws : Applications to structured cell populations[END_REF]). These doubling transmission conditions are handled thanks to discontinuous flux conditions at the corresponding interfaces. Such conditions are tricky to discretize numerically, even on a uniform grid, and require a specific scheme in order to achieve high order of accuracy. We refer the reader to [START_REF] Aymard | A numerical method for transport equations with discontinuous flux functions: Application to mathematical modeling of cell dynamics[END_REF] for a complete study of a Finite Volume (FV) scheme well suited to this purpose. We validate our adaptive strategy on both the localized and distributed mitosis models. The specific application we have in mind is the numerical calibration of terminal follicular development modeled in [START_REF] Echenim | Multi-scale modeling of the follicle selection process in the ovary[END_REF]. Realistic simulations involve solving simultaneously about twenty coupled PDEs in two space dimensions, during about one dozen of cell cycles, with an initial cell population concentrated in the first cell cycle. The weak structure of the coupling is suited to an efficient implementation on parallel architecture, already described in [START_REF] Aymard | Numerical simulation of the selection process of the ovarian follicles[END_REF]. Nevertheless one simulation takes several minutes to run and the extensive simulations required to calibrate the model therefore motivate our search for additional reduction of computational costs by using a MR adaptive finite volume strategy. In this paper we address the problem arising from coupling the MR strategy with the FV method designed in [START_REF] Aymard | A numerical method for transport equations with discontinuous flux functions: Application to mathematical modeling of cell dynamics[END_REF] in the scalar case and with a single processor architecture, leaving the vector case with parallel implementation to future works. Actually, the main difficulty arising from this coupling is already present in the scalar case, where the resulting adaptive scheme (FVMR) unfortunately exhibits unstable numerical noise. The mesh adaptation is tuned by thresholding the details of the MR analysis, which introduces oscillations. In the standard conservative set-up, this spurious noise is bounded by a small parameter, and remains so at future times. In our non conservative situation it is amplified by the discontinuous flux conditions and becomes apparent only after long time simulations -let us say more than ten cell cycles. This corresponds to a proliferation ratio of 2 10 which explains how any numerical error, well within any given tolerance at the beginning of the simulation, can end up having macroscopic consequences on the solution mesh or even its shape (see Figure 14). It is a quite specific set-up, obviously not encountered in all PDEs with discontinuous fluxes. Hopefully they do not all require any particular adaptation of the MR method. This is for instance the case for Multiresolution applied to a traffic flow model [START_REF] Bürger | Fully adaptive multiresolution schemes for strongly degenerate parabolic equations with discontinuous flux[END_REF] with only one (fixed and a priori known) flux discontinuity. If the discontinuity is well captured by the initial graded tree, it will remain that way, which explains why no specific strategy is needed. We will show that the non conservative nature of the equation is indeed the factor responsible for the noise amplification.

The main object of this article is therefore twofold: first to fully document the unstable noise appearance when using the standard MR method on such an equation with doubling conditions. This is done in section 2 where we introduce a simple 1D model and detail the elementary numerical steps, showing how they lead to spurious numerical noise. Secondly, we propose a numerical strategy to stop the noise appearance, and then validate it by extensive numerical tests. This numerical validation, exposed in section 3, uses the FVMR method that we also use in the latest section 4 to solve more realistic biological models. Eventually we show that our method correctly handles both types of modelslocalized or distributed mitosis -with CPU and memory gain between 10 and 20 for an overall relative error remaining below 10 -3 .

Investigation of numerical noise appearance

In this section we investigate the circumstances under which numerical noise can appear. For this we consider a simplified 1D problem with piecewise constant speeds. We compute the approximated solution with a simple first order FV scheme. This allows us to perform one time step of the algorithm by hand and completely describe the interaction of the doubling condition with a two-level MR algorithm. We illustrate this example by numerical simulations performed with a method coupling the 2D FV scheme with a mean value MR adaptive strategy.

Description of the simplified set-up

We study the simple case of a 1D transport model with piecewise constant speeds and discontinuous flux conditions

∂ t φ(t, x) + ∂ x (g(x)φ(t, x)) = 0, for x ∈ R, t > 0 φ(0, x) = φ 0 (x), for x ∈ R, (1)
with N c interfaces numbered from 1 to N c spaced by x s . The speeds are constant in each subdomain

g(x) = g p-1 for (p -1)x s ≤ x ≤ px s , p = 1, . . . , N c + 1.
Transmission conditions are defined by

g p φ(t, px + s ) = k p g p-1 φ(t, (p -1)x - s ) (2)
where k p is equal to 2 to model a doubling flux interface, or to 1 to model a continuous flux interface.

For the sake of simplicity we first choose g p = 1 and k p = 2 for all p = 0, . . . , N c and we apply a doubling condition on the outer boundary

φ(t, 0 + ) = 2φ(t, N c x - s ).
This model is a simple prototype for transport equations modeling cell proliferation with localized mitosis at the end of each cell cycle, which is relevant for the granulosa cells in ovarian follicles [START_REF] Echenim | Multi-scale modeling of the follicle selection process in the ovary[END_REF]. We will also consider an alternative model frequently encountered in the literature (see for instance [START_REF] Perthame | Transport Equations In Biology[END_REF] chapter 3), which consists in distributing the mitosis over the cell cycle duration, with a linear source term

∂ t φ(t, x) + ∂ x (g(x)φ(t, x)) = Bφ(t, x), for x ∈ R, t > 0, φ(0, x) = φ 0 (x), for x ∈ R. (3)
The coefficient B in the source term must be such that the averaged total mass is doubled at the end of a cell cycle, which leads to

exp(BT ) = 2,
where T is the duration of the cycle.

Finite volume discretization on a two-level MR hierarchy

We start by discretizing the solution on two dyadically embedded levels ∆x 0 = 2∆x where we denote the meshes on the finer level by

D 1,k = [k∆x, (k + 1)∆x],
and the meshes on the coarser level by

D 0,k = D 1,2k ∪ D 1,2k+1 = [k∆x 0 , (k + 1)∆x 0 ].
In our double index notation , k, the first number indicates the level, = 0 for the coarser and = 1 for the finer. We further assume that x s = 2K∆x = K∆x 0 , meaning that the p-th doubling interface is located between meshes D 0,pK-1 and D 0,pK on the coarser level, and between meshes D 1,2pK-1 and D 1,2pK on the finer level, as illustrated in Figure 1 for the first interface (p = 1).

3 0,K ∆ ∆ D 0,K-2 D 0,K+1 x x 2 x2 D 1,2K-2 D 1,2K D 1,2K-1 D 1,2K+1 D 0,K-1 D
Figure 1: Two level discretization with doubling interface located between meshes D 0,K-1 and D 0,K on the coarser level, or D 1,2K-1 and D 1,2K on the finer level.

The FV approximation of the solution is defined on the finer level by

φ n 1,k ≈ 1 ∆x D 1,k φ(n∆t, x)dx,
and the solution on the coarser level is defined as

φ n 0,k = 1 2 φ n 1,2k + φ n 1,2k+1 . (4)
Knowing the mean values of the solution on the coarser level, approximated mean values on the finer level can be obtained by local reconstruction. Arbitrary orders of polynomial approximation can be used but we will consider in the sequel only the zeroth order reconstruction, also known as Haar reconstruction in the wavelet context

φn 1,2k = φn 1,2k+1 = φ n 0,k , (5)
and the second order reconstruction

φn 1,2k = 1 8 φ n 0,k-1 + φ n 0,k - 1 8 φ n 0,k+1 , φn 1,2k+1 = - 1 8 φ n 0,k-1 + φ n 0,k + 1 8 φ n 0,k+1 . (6)
We use a FV scheme to compute the numerical solution as

φ n+1 j = φ n j - ∆t ∆x j F n j+1/2 -F n j-1/2 , (7)
where F n j+1/2 is the numerical flux at time t n through the edge between meshes D j and D j+1 , and ∆x j the size of mesh D j . A simple upwind first order scheme on a uniform grid consists, in this case with unit speed, in taking

F n j+1/2 = φ n j if j = 0, pK -1 and 1, 2pK -1, 2φ n j otherwise, ( 8
)
where the factor 2 in bold takes into account the doubling flux condition (2) when k p = 2, for p = 1, . . . , N c.

We can therefore compute the numerical solution on the finer level as

φ n+1 1,k = (1 -λ)φ n 1,k + λφ n 1,k-1 if k = 2pK, (1 -λ)φ n 1,k + 2λφ n 1,k-1 if k = 2pK, for p = 1, . . . , N c, (9)
with λ = ∆t ∆x , which must satisfy 0 < λ ≤ 1 in order to ensure the stability of the solution. If the solution is smooth enough -with a criterion that we will precise below -we can compute it instead on the coarser level, using λ/2 instead of λ

φ n+1 0,k = (1 -λ 2 )φ n 0,k + λ 2 φ n 0,k-1 if k = pK, (1 -λ 2 )φ n 0,k + λ 2 2φ n 0,k-1 if k = pK, for p = 1, . . . , N c.
Note that the numerical solution becomes locally small but nonzero, even when starting from a piecewise constant solution: due to the numerical diffusion of the numerical scheme, the solution rapidly becomes nonzero upstream of the initial condition. For instance for the simple scheme [START_REF] Coquel | Convergence of time-space adaptive algorithms for nonlinear conservation laws[END_REF] we have

φ n 1,k = (1 -λ)φ n-1 1,k + λφ n-1 1,k-1 , = (1 -λ) 2 φ n-2 1,k + 2(1 -λ)λφ n-2 1,k-1 + λ 2 φ n-1 1,k-2 , . . . = . . . = (1 -λ) n φ 0 1,k + . . . + λ n φ 0 1,k-n . (10)
In the general case with an arbitrary number of levels and arbitrary high dimensions, one time step of the FVMR scheme consists in a loop over levels starting from the finest. For each level, a first loop computes the numerical fluxes between meshes active at the current level using [START_REF] Coquel | Local time stepping applied to implicit-explicit methods for hyperbolic systems[END_REF], possibly requiring local reconstruction if neighboring meshes entering the flux estimation are coarser than the current level. Then a second loop updates the solution on the meshes of the current level using [START_REF] Coquel | Multiresolution and Adaptive Methods for Convection-Dominated Problems[END_REF].

Having in mind this algorithm is necessary to correctly handle the cases of neighboring meshes of different sizes

• if the left hand side neighbor of mesh D 1,2k is the coarser mesh D 0,k-1 , then the solution φn 1,2k-1 at time n on mesh D 1,2k-1 must be reconstructed, using [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] or [START_REF] Cohen | Adaptive multiresolution for finite volume solutions of gas dynamics[END_REF], in order to compute F n (1,2k)-1/2

• if the left hand side neighbor of mesh D 0,k is the finer mesh D 1,2k-1 , then

F n (0,k)-1/2 = F n (1,2k-1)+1/2 = φ n 1,2k-1
has already been computed during the loop on the finer level

• if the right hand side neighbor of mesh D 1,2k+1 is the coarser mesh D 0,k+1 , no reconstruction is required to compute F n (1,2k+1)+1/2 for the simple first order scheme [START_REF] Coquel | Local time stepping applied to implicit-explicit methods for hyperbolic systems[END_REF]. In the case of a higher order scheme, φn

1,2k+2
would have to be reconstructed

• if the right hand side neighbor of mesh D 0,k is the finer mesh D 1,2k+2 , then

F n (0,k)+1/2 = F n (1,2k+2)-1/2 = φn 1,2k+1
has already been computed during the loop on the finer level, using reconstructed value φn

1,2k+1

The choice between the two levels of discretization on a given mesh D 0,k in the domain is made by comparing the solution on its subdivisions D 1,2k ∪ D 1,2k+1 with the solution reconstructed on the same mesh using the values on the coarser grid with formulas (5) or [START_REF] Cohen | Adaptive multiresolution for finite volume solutions of gas dynamics[END_REF]. The local smoothness of the solution is indeed measured by the detail, defined by

d n 1,k = φ n 1,2k -φn 1,2k (11)
• In the Haar case, a null detail corresponds to a locally constant function

d n 1,k = φ n 1,2k -φ n 0,k = 1 2 (φ n 1,2k -φ n 1,2k+1 ). (12)
• In the quadratic basis, a null detail corresponds to a locally quadratic function

d n 1,k = φ n 1,2k - 1 8 φ n 0,k-1 + φ n 0,k - 1 8 φ n 0,k+1 = - 1 16 (φ n 1,2k-2 + φ n 1,2k-1 ) + 1 2 (φ n 1,2k -φ n 1,2k+1 ) + 1 16 (φ n 1,2k+2 + φ n 1,2k+3 ). (13)
We refer for instance to [START_REF] Harten | Multiresolution algorithms for the numerical solutions of hyperbolic conservation laws[END_REF] for the approximation theory results establishing the link between d n 1,k and the local smoothness of the solution.

In order to monitor the adaptive FV scheme, the absolute value of the details are tested at each time step against a level dependent threshold ε ε = ε2 -L+1 , for = 1, . . . , L -1.

(

) 14
In the numerical simulations, we refer to the MR threshold to indicate the threshold value on the finest MR level ε = ε L-1 . The set of indices of the significant details at time t n is designed as

Γ n ε = {( , k), ||d ,k || ≤ ε},
and enlarged into a set Γn+1 ε according a heuristic rule [START_REF] Harten | Multiresolution algorithms for the numerical solutions of hyperbolic conservation laws[END_REF], accounting for singularity displacement (and formation in the nonlinear case) between time t n and t n+1 . Here in the two level case, if |d n 1,k | ≥ ε 1 , then not only D 0,k but also the two neighbors D 0,k-1 and D 0,k+1 should be subdivided to compute the solution at the next time step. The whole predictive procedure of adaptive grid refinement is summarized in the general case with L MR levels in Algorithm 1. The solution (φ n ,k ) ( ,k)∈ §n+1 is reconstructed on the adaptive grid at any given time t n with Algorithm 2. It is then updated to time t n+1 with a FV scheme. Theoretically and practically -although it implies to renounce the MR computational benefits-it is also possible to compute the reconstructed solution φ n ε on the finest grid using Algorithm 2 and setting all d n l,k = 0 for (l, k) ∈ Γn+1 0 \ Γn+1 ε . In the case of a scalar hyperbolic equation numerically discretized with a Lipschitz flux, the L 1 norm error a n between φ n ε and the FV solution on the finest level φ n unif has been shown in [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] to be bounded

a n := ||φ n unif -φ n ε || ≤ C(t n )ε. (15)
The threshold ε should be chosen so that a n is of the same order as the estimate of the discretization error e n

e n := ||φ n unif -φ n exact || ≤ C(t n )∆x p , ( 16
)
where p is the numerical order of the reference FV scheme on the uniform grid. Of course in realistic set-ups the exact solution is unavailable and one also wishes to avoid computing the FV reference solution on the fine grid. Estimates ( 15) and ( 16) can be used to calibrate the threshold value on a smaller scale domain in time and space, or for a simplified model. In our case, we will see in the numerical tests that choosing ε between 10 -2 and 10 -3 supplies a good compromise between accuracy and efficiency.

Explication of the numerical noise appearance

In this paragraph we describe the mechanism for numerical noise appearance due to the coupling of Multiresolution with doubling fluxes.

Algorithm 1 Prediction of the set Γ n+1 ε of significant details at time t n+1 from the details at time t n . General case with L MR levels numbered from the coarsest 0 to the finest L -1 Require: The solution is known by its mean values on the coarsest level and the details {d n ,k } ( ,k)∈Γ n ε Prediction:

for level = L -1 1 do for k = 0 N 0 2 -1 do if ( , k) ∈ Γ n ε and |d n ,k | ≥ ε then ( , k + i) ∈ Γ n+1 ε for |i| ≤ 1 if < L -
for level = L -1 1 do for k = 0 N 0 2 -1 do if ( , k) ∈ Γ n+1 ε then for |i| ≤ 1 do ( -1, k/2 + i) ∈ Γ n+1 ε end for end if end for for k = 0 N 0 2 -1 do if ( , k) ∈ Γ n+1 ε \Γ n ε then d n ,k = 0 end if end for end for Suppose that |d n K-i | < ε, for i ≥ 0 and |d n K+1 | ≥ ε.
Due to the preceding rule, the mesh D 0,K = D 1,2K ∪ D 1,2K+1 is subdivided, to account for possible displacement of the singularity (see Algorithm 1). The mesh D 0,K-1 is kept at the coarser level, but the value of the solution on its right subdivision D 1,2K-1 , which enters the FV scheme to advance the solution on D 1,2K , is available using the prediction operator [START_REF] Cohen | Adaptive multiresolution for finite volume solutions of gas dynamics[END_REF] in the quadratic case.

Let us compute the solution at time n + 1 in the neighborhood of the doubling interface (see Figure 1)

φ n+1 0,K-1 = φ n 0,K-1 - λ 2 φn 1,2K-1 -φ n 0,K-2 , φ n+1
1,2K = (1 -λ)φ n 1,2K + 2λ φn 1,2K-1 , φ n+1 1,2K+1 = (1 -λ)φ n 1,2K+1 + λφ n 1,2K , φ n+1 1,2K+2 = (1 -λ)φ n 1,2K+2 + λφ n 1,2K+1 , φ n+1 1,2K+3 = (1 -λ)φ n 1,2K+3 + λφ n 1,2K+2 .
Algorithm 2 Construction of the solution on the adaptive grid S n+1 from the set of significant details Γ n+1 ε Require: The solution is known by its mean values on the coarsest grid, and the details for

(k, ) ∈ Γ n+1 ε for level = 1 L -1 do for k = 0 N 0 2 -1 do if ( , k) ∈ Γ n+1 ε then {D ,2k , D ,2k+1 } ∈ S n+1 and D -1,k / ∈ S n+1 compute φ n ,2k = φn ,2k + d n ,k and φ n ,2k+1 = φn ,2k+1 -d n ,k

end if end for end for

We should now check the smoothness of the solution after this update. Using [START_REF] Louvet | Summer school on multiresolution and adaptive mesh refinement methods[END_REF] and ( 4) we can compute d n+1 K as

d n+1 1,K = φ n+1 1,2K -φn+1 1,2K = φ n+1 1,2K - 1 8 φ n+1 0,K-1 + φ n+1 0,K - 1 8 φ n+1 0,K+1 = - 1 8 φ n+1 0,K-1 + 1 2 (φ n+1 1,2K -φ n+1 1,2K+1 ) + 1 16 (φ n+1 1,2K+2 + φ n+1 1,2K+3 ).
We plug in the expressions of the solution at time n + 1 and obtain

d n+1 K = - λ 16 φ n 0,K-2 + λ 2 (2 φn 1,2K-1 -φ n 1,2K ) + λ 16 (φ n 1,2K+1 + φ n 1,2K+2 ) - 1 8 φ n 0,K-1 + λ 16 φn 1,2K-1 + 1 -λ 2 (φ n 1,2K -φ n 1,2K+1 ) + 1 -λ 16 (φ n 1,2K+2 + φ n 1,2K+3 ).
The effect of the doubling flux on the detail is an additional term

d n+1 1,K = dn+1 1,K + λ 2 φn 1,2K-1 , (17)
compared to the detail dn+1 K if no doubling condition is imposed. This means that if the solution is smooth at time t n , but the solution is large enough before the doubling interface, then the detail may become significant due to the doubling.

We now explain what happens next, if the detail d n+1

1,K which was negligible at time t n becomes larger than the MR threshold ε at time n + 1. In the quadratic case, as described in Algorithm 1, the mesh D 0,K-1 just upstream of the doubling interface should be predictively refined using ( 6)

φn+1 1,2K-2 = 1 8 φ n+1 0,K-2 + φ n+1 0,K-1 -1 8 φ n+1 0,K , φn+1 1,2K-1 = -1 8 φ n+1 0,K-2 + φ n+1 0,K-1 + 1 8 φ n+1 0,K ,
where

φ n+1 0,K = 1 2 (φ n+1 1,2K + φ n+1 1,2K+1 ), = 1 -λ 2 (φ n 1,2K + φ n 1,2K+1 ) + λ 2 (2 φn 1,2K-1 + φ n 1,2K ). (18)
Here again the difference between the values of the solution on the subdivisions, which should be small since the solution is only predictively refined, is enhanced by the value φn 1,2K-1 upstream of the interface which contributes twice instead of once. We see that in a single time step the doubling condition coupled with the MR, has introduced an oscillation of the order of φn 1,2K-1 .

Suppose we start from an initial condition localized in the first slab before the first doubling interface, for instance φ 0 k = δ 0,k . Due to the numerical diffusion effects described earlier (see [START_REF] Echenim | Multiscale modeling of follicular ovulation as a reachability problem[END_REF]), after M time steps (M depending on ε and λ), the solution φ M k is nonzero and smooth for k = 0, . . . , M . The MR scheme will therefore coarsen the grid in this region. As the solution is transported over several slabs, its amplitude is doubled at each interface, and depending on its local amplitude we can fall in the situation depicted by [START_REF] Sakaue-Sawano | Visualizing spatiotemporal dynamics of multicellular cell-cycle progression[END_REF]. Furthermore, another factor of noise appearance is the dependence of the MR threshold level on the discretization level according to the rule [START_REF] Michel | Multiscale modeling of follicular ovulation as a mass and maturity dynamical system[END_REF].

We see from [START_REF] Sakaue-Sawano | Visualizing spatiotemporal dynamics of multicellular cell-cycle progression[END_REF] that the decision to refine the mesh after crossing a doubling interface depends not only on the smoothness but also on the local amplitude. Thanks to ( 14), the same (small) amplitude discretized on a coarse level will be more likely to trigger the refinement, and might even trigger the local refinement on two consecutive levels (see Algorithm 1). We propose an example for such a situation in Appendix B. We describe a situation with three levels of discretization, where the solution has sensibly the same small amplitude but is refined partly on the intermediate level, and partly on the coarsest one. We show that refinement rules together with rules ( 17) and ( 14) lead to coarsen the solution wherever it was previously on the intermediate level, and to refine it on the finest level on some locations where it was before on the coarsest one. At this point we also stress out that the difference between the two wavelets ( 12) and ( 13) is quite important, specially if we solve the 1D test case presented above with the 2D FV code. To illustrate this point we perform a simulation on a domain consisting of twelve subdomains of unit width and height with vertical doubling interfaces. The (horizontal) speed is equal to 1. The initial condition is piecewise constant φ(t = 0, x, y) = 10 for 0 ≤ x < 1 and 0.2 ≤ y < 0.3 0 elsewhere

We first use a two level MR hierarchy, with five meshes per subdomain in the x direction on the coarser level. that due to the detail amplification ( 17) the MR solutions with both Haar and quadratic reconstructions exhibit more refined meshes upstream than downstream of the support of the solution. Figure 3 shows the solution with a color code from 0 to 10 corresponding to the initial condition; the spurious noise is visible at this scale. Note that in the Haar case, the upstream amplification ( 18) does not happen. In the case of the quadratic reconstruction, solving the 1D transport equation ( 1) with the 2D MR scheme also generates some nonzero solution in the orthogonal direction, which gets progressively amplified by the doubling interfaces. The difference between the panels corresponding to the Haar (panel a)) and quadratic predictions (panel b)) highlights the appearance of noise due to the quadratic prediction operator in 2D: the solution becomes nonzero in the transversal direction above and below the support of the initial condition. The quadratic prediction performs a better compression of the solution than the Haar one, therefore the mesh is coarsened faster away from the support of the solution. However in both cases the details are later amplified by [START_REF] Sakaue-Sawano | Visualizing spatiotemporal dynamics of multicellular cell-cycle progression[END_REF] and, in the quadratic case, also by [START_REF] Shang | Cauchy problem for multiscale conservation laws : Applications to structured cell populations[END_REF], which explains the resurgence of refined meshes upstream of the support of the solution (panels b)). In this simple set-up isolated cells appear on the finest level only in the quadratic case, but in more complicated situations the same phenomenon also occurs in the case of the Haar prediction.

Adaptive doubling strategy

To prevent the formation of this "numerical noise" several strategies can be implemented :

• Forcibly refining the solution on the finest discretization grid in the vicinity of interfaces, without considering its local smoothness. This may become quite greedy and altogether compensate the benefits of the MR strategy.

• Adapting the MR threshold to the average amplitude of the solution. This should be done carefully to avoid rapid changes in the threshold and subsequently rapid coarsening of the solution. We will document later that the optimal strategy is to adapt the threshold to the total mass of the solution

ε t, = ε M (t) (19)
where the level dependent ε is defined by ( 14), with ε calibrated for a unit mass solution at initial time.

• Modifying the transmission condition (2) with an adaptive strategy consisting in applying the doubling only if the solution is above a given threshold δ

g p φ(t, px + s ) = 2g p-1 φ((p -1)x - s ) if |φ(t, (p -1)x - s )| ≥ δ, g p-1 φ((p -1)x - s ) otherwise. ( 20
)
In practice δ should be chosen close to ε, and we will study the relationship more closely in paragraph 3.3.

We also assume that the support of the initial condition is entirely located in the first slab before the first doubling interface. We have actually exemplified in Appendix B that a small and smooth solution lying on both sides of an interface, could right away generate spurious oscillations. Restricting the initial support ensures that no additional oscillations will be created by coupling the MR with the doubling condition. The bottom panels of Figures 2 and3 display the effect of the adaptive doubling strategy in the case of second order reconstruction, setting the parameter δ equal to ε = 10 -2 . The numerical noise has completely disappeared.

Formalization and generalization of the doubling strategy

The result ( 17) also applies to Haar reconstruction [START_REF] Harten | Multiresolution algorithms for the numerical solutions of hyperbolic conservation laws[END_REF]. This means that the appearance of the spurious noise is due not only to the discontinuities in the solution induced by the localized doubling, but also to the enhancement of the thresholding error. In particular, even if we modified the prediction scheme in the vicinity of interfaces in order to account for the doubling of the solution, this would not stop the increase of noise.

To formalize this remark, we generalize the MR scheme introduced in Paragraph 2.2 to an arbitrary number of levels L. We introduce the linear prediction operator P -1 such that formula ( 6) can be expressed generically as

φn = P -1 φ n -1 ,
and the reconstruction formulas in Algorithm 2 as

φ n = P -1 φ n -1 + d . (21)
Operator P -1 is a 2 N 0 × 2 -1 N 0 matrix, that is zero everywhere except at three consecutive coefficients per row

(P -1 ) 2k,k-1 = 1/8, (P -1 ) 2k,k = 1, (P -1 ) 2k,k+1 = -1/8, (P -1 ) 2k+1,k-1 = -1/8, (P -1 ) 2k+1,k = 1, (P -1 ) 2k-1,k+1 = 1/8,
with periodic boundary conditions

(P -1 ) 0,2 -1 N0-1 = 1/8, (P -1 ) 2 N0-2,0 = -1/8, (P -1 ) 1,2 -1 N0-1 = -1/8, (P -1 ) 2 N0-1,0 = 1/8. (22)
Note that other types of boundary conditions can be designed but periodic ones are coherent with the boundary conditions used in the application to the biological model. Formula (21) can be applied recursively, leading to

φ n L-1 = L-1 =1 P -1 φ n 0 + L-1 =1 L-1 j= +1 P j-1 j d n . (23)
The adaptive solution after thresholding with the MR threshold ε can be expressed as

φ n,ε L-1 = L-1 =1 P -1 φ n 0 + L-1 =1 L-1 j= +1 P j-1 j d n ,ε , (24)
where d n ,ε is the vector of thresholded details on level

(d n ,ε ) ,k = d n ,k if ||d n ,k || ≥ ε , 0 otherwise.
From ( 23) and (24) it is straightforward to obtain an O(ε) estimate of the error between the unthresholded solution φ n L-1 and the thresholded one φ n,ε L-1 . In [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF], a strategy was introduced to enlarge the tree of significant details (with respect to the threshold ε), ensuring that the adaptive solution on the corresponding mesh satisfies the same error estimate at the next time step time n + 1. We also recall a useful result from [START_REF] Coquel | Convergence of time-space adaptive algorithms for nonlinear conservation laws[END_REF] showing that the reconstructed solution on any intermediate level is vitiated by an O(ε) perturbation due to MR thresholding and reconstruction, which we can summarize as

φ n,ε ,k = φ n ,k + O(ε). (25)
We apply the adaptive doubling strategy (20) along with the elementary FV scheme [START_REF] Coquel | Convergence of time-space adaptive algorithms for nonlinear conservation laws[END_REF] generalized to an L level MR hierarchy with

λ = 2 -L+1 λ φ n+1,ε ,k = (1 -λ )φ n,ε ,k + λ φ n,ε ,k-1 if k = pK or |φ n,ε ,k-1 | ≤ δ (1 -λ )φ n,ε ,k + λ 2φ n,ε ,k-1 if k = pK and |φ n,ε ,k-1 | > δ,
with p = 1, . . . , N c . We see that we have a conditional maximum principle in the sense that

if φ n,ε ,k ≤ δ and φ n,ε ,k-1 ≤ δ then φ n+1,ε ,k ≤ δ. (26)
Relations ( 25) and ( 26) together show that parameter δ should be correctly calibrated in relation with the MR threshold ε, to prevent the MR perturbation in (25) from being amplified and becoming the unwanted spurious noise.

Numerical validation of the adaptive thresholds strategies

In this paragraph we perform numerical tests to validate the adaptive strategies ( 19) and (20). We then study the relationship between the two thresholds δ and ε.

Effect of the adaptive doubling strategy on the solution

To justify the doubling strategy designed to tackle the spurious noise introduced by the MR we show the influence of applying (20) instead of (2), independently of the MR, and even of any numerical scheme, by computing the exact solution of the 1D problem [START_REF] Aymard | Numerical simulation of the selection process of the ovarian follicles[END_REF]. As shown in [START_REF] Aymard | A numerical method for transport equations with discontinuous flux functions: Application to mathematical modeling of cell dynamics[END_REF], the exact solution for a single interface located at x = x s with condition (2) is

φ 1 (t, x) =      φ 0 (x -g 0 t) for x < x s , φ1 (t s ) for (x s -g 1 t) ≤ (x -g 1 t) < x s , φ 0 (x -g 1 t) for (x -g 1 t) ≥ x s ,
where φ1 is the trace of the solution on the right of the interface, defined by

φ1 (t) = 1 g 1 k 1 (g 0 φ 0 (x s -g 0 t)), (27)
and t s is the delay after which the effect of interface x s is felt at position x and time t defined by

t s = t - (x -x s ) g 1 . ( 28
)
We now show how the exact solution for a sequence of N c interfaces can be defined recursively. Let us denote by φ p-1 (t, x) the solution for a sequence of p -1 interfaces, with p > 1 and add a p th interface at px s with a coefficient k p for the transmission condition and a speed g p on the right hand side (see Figure 4). For x < px s the effect of the p th interface is not felt and the solution is φ p-1 (t, x). For x ≥ px s , depending on time t, the effect of the interface is felt or not. If (x -px s ) < g p t, events crossing the p th interface are multiplied by the k p coefficient and transported at speed g p . The effects of these events are felt at position x after a delay (x -px s )/g p . If (x -px s ) > g p t, the effect of the p th interface is not felt yet and the initial condition defined on x > px s is transported at speed g p . We eventually have

φ p (x, t) =      φ p-1 (x, t) for x < px s , φp (t p s ) for (px s -g p t) ≤ (x -g p t) < px s , φ 0 (x -g p t) for x -g p t ≥ px s , (29)
where φp is the trace of the solution on the right of the p th interface, defined by

φp (t) =      1 g p k p g p-1 φ p-1 (px s , t) if p > 1 1 g p k p g p-1 φ 0 (x s -g p-1 t) otherwise ( 30
)
and t p s is defined by

t p s = t - (x -px s ) g p .
We mimic set-up (1) by solving equation (29) numerically using recursiveness, with g p = 1 for all p and an initial condition φ 0 defined by

φ 0 (x) = exp -5(x -0.5) 2 . ( 31
)
We also mimic strategy (20), that we plan to apply to our real problem, and use the doubling condition in (30) only if the solution is larger than the doubling threshold δ. At each time step, we compute the total mass by numerically integrating the solution over the whole x abscissa range [0, 20].

Err(δ) = n M δ,n -M n n |M n | , ( 32
)
where φ δ (respectively φ) is computed on the discrete time-space grid ∆t × ∆x with (29) and with (respectively without) the adaptive doubling strategy (20). The error is displayed as a function of the threshold parameter δ, in logarithmic axes, which highlights its linear dependence onto the threshold. In Figure 6 we also display the error behavior for the same simulation in the case when the parameter δ changes with time, which would happen if the MR threshold ε followed the time-dependent rule [START_REF] Tomura | Contrasting quiescent G0 phase with mitotic cell cycling in the mouse immune system[END_REF]. We choose in fact

δ t = δM (t), (33)
where M (t) is the mass of the solution integrated over the whole spatial range, at time t

M (t) = N xs 0 φ(t, x)dx.
These numerical experiments show that whichever the selected thresholding strategy for the MR method, the influence of the adaptive doubling strategy (20) is well controlled by parameter δ. Except for the unreasonably large value of δ = 0.1, the mass behavior computed with the adaptive doubling strategy is undistinguishable from the exact solution.

Both error curves in Figure 6 exhibit a better behavior than the asymptotic O(δ).

In Figure 7 we check that the numerical scheme FVMR described in the previous section correctly reproduces this macroscopic mass behavior. We run a 2D simulation of the 1D simple model ( 1) and we compare it to case (3) where the doubling is uniformly distributed over the cycle, with a linear birth term amounting to the same overall doubling.

In Figure 8 snapshots of solutions at t = 7.48 and t = 8 are displayed, allowing one to compare at the microscopic level the effect of either applying a doubling condition locally at the end of interface or smoothly distributing proliferation over all the cycle. At both times the solutions sensibly have the same mass, we have selected the snapshot at times when the red and dashed green curves intersect in Figure 7. The top panel, at t = 7.48, corresponds to a time when the support of the solution is equally distributed on one side and the other of an interface between cell cycles and the induced discontinuity in the solution is clearly visible on the left panel. The bottom panel, at t = 8 corresponds to a time when the support of the solution is well within the same cell cycle, and both solutions do not substantially differ.

Time adaptation of the MR threshold ε t

To study the influence of the time adaptation of the MR threshold ε t [START_REF] Tomura | Contrasting quiescent G0 phase with mitotic cell cycling in the mouse immune system[END_REF], we perform different simulations of simple test case (1) with smooth initial condition (31), the third order FV scheme and a 4 level MR hierarchy with the quadratic reconstruction. We test by turns, the time adaptation of the MR threshold [START_REF] Tomura | Contrasting quiescent G0 phase with mitotic cell cycling in the mouse immune system[END_REF], and another adaptation of the MR threshold based on the maximum value of the solution

ε t, = ε max x,y |φ(t, x, y)|. (34)
In each case we perform the simulation with or without the adaptive doubling strategy (20). The six possible cases are referred in the graphs with the following legends Unconditional doubling (2) Adaptive doubling (20)

Constant MR threshold δ = 0 × δ = ε + Adaptive MR threshold (19) δ = 0, ε t = εM (t) δ = ε t = εM (t) * Adaptive MR threshold (34) δ = 0, ε t = εSup|Φ|(t) • δ = ε t = εSup|Φ|(t)
Table 1: Legend codes for Figures 9 and15 The L 1 -norm of the relative error at a given time t n = n∆t is defined as

||E|| L 1 = i,j φ ε,n L-1;i,j -φ n L-1;i,j i,j φ n L-1;i,j , (35)
where φ n L-1;i,j is the FV solution at time t n = n∆t on the finest MR level L -1 of discretization, with index i in the x direction and j in the y direction, and φ ε,n L-1;i,j is the adaptive scheme solution, reconstructed at time t n on the finest level L -1.

This error is computed at five different times and displayed in Figure 9 against the MR threshold level ε on the left panels, the CPU gain on the center panels, and the memory gain on the right panels. The CPU gain is the ratio of the CPU time for the simulation on the uniform finest grid to the CPU time for the FVMR simulation. The memory gain is the ratio of the number of meshes in the uniform finest grid to the maximum number of meshes in the adaptive grid during the FVMR simulation.

The left panels show that, whichever the thresholding strategy, the expected asymptotic behavior in O(ε) is satisfied (see equation ( 15)). Nevertheless, as time increases, the normalization factors in the adaptive strategies [START_REF] Tomura | Contrasting quiescent G0 phase with mitotic cell cycling in the mouse immune system[END_REF] and (34) increase, resulting in more and more distinct curves. In practice, this type of curve is computed for a reduced version of the full scale simulation in order to select a threshold value ε ensuring a correct balance between the MR (15) and discretization ( 16) errors. In our case we will use the optimal values of the thresholds δ and ε calibrated from this simple 1D set-up in future 2D and vector simulations of the biological model. In other situations, one could calibrate the MR threshold in 2D tests before addressing the 3D cases where the reference solution is beyond reach. We have represented the error in L 1 norm in space at various times during the simulation, which gives more detailed information than the global L 1 error over the whole timespan, and allows us in particular to see the variation in the error as the time adaptive threshold increases. For a constant threshold (red and green lines) the error remains almost unchanged with time. This is not a surprise since we have already observed in [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] that the theoretical estimate (15) predicting a linear growth in time was very pessimistic. Also quite expectedly, the difference between the error curves increases with time in the case of time adaptive thresholding with increasing value of ε. Note that the MR error at a given time has built up progressively and does not entirely depend on the current value of ε(t): at time t = 10 there is a factor of order 10 between the red and dark blue curves, which correspond respectively to the constant and mass dependent threshold (Eq. ( 19)and (33)), although at that time the mass, and therefore ε(t), has already been multiplied by 10 3 . At the end of the simulation (at t = 20) the error ratio between both curves is still about 10, although ε(t) has been multiplied by 10 6 .

The center and right panels allow one to assess the performance of the thresholding strategies with respect to computational costs, in terms of CPU time as well as memory requirement. For early times all strategies seem equivalent, but after ten cycles, the simulations corresponding to constant MR thresholding start having worse performances than those using either of adaptive strategies ( 19) and (34), which remain basically undistinguishable. The simulation with constant MR thresholding and no doubling adaptation, displayed with green ×, has the worst performances, because at large times, spurious noise has started to develop, requiring fine mesh discretization and therefore, more memory and more CPU time.

We now validate the adaptive doubling strategy (20) by extensive numerical tests. We compare its efficiency when it is coupled with one or another method of MR thresholding. We also estimate the optimal value for the threshold parameter δ.

Link between the MR and the doubling thresholds

In this paragraph we study the relation between the MR threshold parameter ε and the threshold δ in the adaptive doubling interface (20). Still for 1D test case (1) with the piecewise constant initial condition, we run extensive tests for different values of both thresholds. We also study the influence of the number of levels in the MR hierarchy, the precision of the FV algorithm, and the effect of normalizing the MR threshold by the overall mass of the solution, as it increases with time (see Eq. ( 19)).

Figure 10 displays the L 1 norm of error at final time t = 20, between the adaptive solution and the solution computed on the uniform finest grid. The adaptive solution is computed with different values of the MR threshold ε, fixed once for all times, and for values of the threshold δ varying between ε/100 and 10ε. For two, three and four levels in the MR hierarchy we display the results in two different fashions : The standard error visualization is shown in the left panels, where the interesting feature is the linear trend of the error with ε which is preserved by the adaptive doubling strategy. However normalizing by ε both the error on the vertical axis and δ on the horizontal axis highlights the 35) at different times, with respect to the MR threshold ε (left panel), the CPU gain (center panel) and the memory gain (right panel). Effect of adaptive MR threshold ε t and doubling threshold δ in the case of the simple 1D test case, 4 MR levels, third order FV scheme. The line labelled O(ε) in the left panel graphs indicates the theoretical asymptotic behavior of the error as a function of ε, which is a straight line of slope one in logarithmic axes. The complete legend code is described in Table 1. correlation between ε and δ since the minimum of the error is observed for a similar ratio range for each set of curves.

Figure 10 displays results computed with the first order upwind FV scheme [START_REF] Coquel | Multiresolution and Adaptive Methods for Convection-Dominated Problems[END_REF][START_REF] Coquel | Local time stepping applied to implicit-explicit methods for hyperbolic systems[END_REF]. The value of δ which minimizes the error lies between 0.2ε and 2ε and increases with both ε and the number of levels. We also made the same test for four levels and the third order FV scheme. The results are displayed in the top panels of Figure 11; the value of δ minimizing the error is shifted from about ε in the first order case to about 2ε.

We then couple the adaptive doubling strategy (20) with the adaptation of the MR threshold to the increasing size of the solution. The same error study as previously is performed, for four MR levels and the third order FV scheme. The results, displayed in the bottom panels of Figure 11, still exhibit that there is an optimal value for δ. It seems to be roughly around 1.5ε.

After these extensive numerical tests, we are confident that coupling both adaptive thresholds [START_REF] Tomura | Contrasting quiescent G0 phase with mitotic cell cycling in the mouse immune system[END_REF] and ( 20) enables the FVMR strategy to robustly handle doubling flux conditions of type (2). We will adopt in the sequel the value δ = ε which is close to the optimal value in all cases.

4 Application to a biological model

Model and numerical method

In this section we apply the strategy presented above to the numerical simulation of a biological model dedicated to the cell dynamics in ovarian follicular development introduced in Echenim et al. [START_REF] Echenim | Multiscale modeling of follicular ovulation as a reachability problem[END_REF] and references therein. The model is multiscale to account for the selection process of ovulatory follicles which involves the cellular, follicular and ovarian levels. At the microscopic level, the granulosa cell population is structured according to the cell age (position within the cell cycle) and the cell maturity (level of sensitivity towards hormonal control). In each ovarian follicle, the granulosa cell population is described by a density function whose changes are ruled by conservation laws. The multiscale structure arises from the formulation of a hierarchical control operating on the aging and maturation velocities as well on the source terms of the conservation law; this control depends on the first moment in the maturity variable of the density. Diagrams summarizing the main features of the model are depicted in Figure 12. The functional space in cell age × cell maturity is subdivided into subdomains corresponding to different cell states: proliferation (Ω p 1 ∪ Ω p 2 ) or differentiation (Ω 3 ), and also with respect to the cell sensitivity towards the hormonal control exerted by the pituitary gland (sensitive phase Ω p 1 ∪ Ω 3 , insensitive phase Ω p 2 ). The mitosis is the endpoint of the cell cycle which is completed when the two daughter cells are separated from each other. In our set-up the mitosis happens at values x s = p for p = 1, . . . , N c -1, on the internal boundaries between Ω p 2 and Ω p+1 1 , and between Ω Nc 2 and Ω 1 1 when the density reaches the right side of the domain due to the periodic boundary condition at this point. On the intermediate internal boundaries between Ω p 1 and Ω p 2 , at values x s = p + 1/2 for p = 1, . . . , N c , the flux condition is continuous. The vertical component of the velocity is positive in Ω p 1 and consequently all cells eventually enter Ω 3 . In the differentiation phase Ω 3 , cells do not proliferate any more and cannot overcome an asymptotical maturity level, corresponding to a null vertical velocity (see [START_REF] Michel | Multiscale modeling of follicular ovulation as a mass and maturity dynamical system[END_REF]). The position of the zero of h(γ, u) (see Eq (40)) varies with the control term u but the properties of function h ensure that it remains in Ω 3 . We can thus apply periodic boundary conditions to the density -which are easy to implement in our adaptive mesh context-on boundaries y = 0 and y = 1. We refer the interested reader to [START_REF] Aymard | A numerical method for transport equations with discontinuous flux functions: Application to mathematical modeling of cell dynamics[END_REF] for a detailed exposition of the motivation underlying the choice of a localized transmission condition of type (2), and to [START_REF] Tomura | Contrasting quiescent G0 phase with mitotic cell cycling in the mouse immune system[END_REF][START_REF] Sakaue-Sawano | Visualizing spatiotemporal dynamics of multicellular cell-cycle progression[END_REF] for biological evidence of localized mitotic phenomenon). We also refer the reader to Appendix A for the description of the model in the single follicle case and to [START_REF] Aymard | A numerical method for transport equations with discontinuous flux functions: Application to mathematical modeling of cell dynamics[END_REF] for a more general presentation, including a description of the third order numerical scheme designed to handle the discontinuous flux conditions. Here we solely address the coupling of the MR with the FV scheme in the context of a localized or distributed mitosis. Actually we even neglect at this stage the macroscopic scale where the different follicles are coupled through the interaction between the ovaries and the hypothalamus/pituitary complex, formerly exposed in [START_REF] Aymard | Numerical simulation of the selection process of the ovarian follicles[END_REF],

||E|| L1 /ε δ/ε ε=10 -9 ε=10 -8 ε=10 -7 ε=10 -6 ε=10 -5 ε=10 -4 ε=10 -3 ε=10 -2 ε=10 -1
Two MR levels, first order FV scheme

||E|| L1 /ε δ/ε ε=10 -9 ε=10 -8 ε=10 -7 ε=10 -6 ε=10 -5 ε=10 -4 ε=10 -3 ε=10 -2 ε=10 -1
Three MR levels, first order FV scheme

||E|| L1 /ε δ/ε ε=10 -9 ε=10 -8 ε=10 -7 ε=10 -6 ε=10 -5 ε=10 -4 ε=10 -3 ε=10 -2 ε=10 -1
Four MR levels, first order FV scheme

||E|| L1 /ε δ/ε ε=10 -9 ε=10 -8 ε=10 -7 ε=10 -6 ε=10 -5 ε=10 -4 ε=10 -3 ε=10 -2 ε=10 -1
Constant MR thresholding ε

||E|| L1 /ε δ/ε ε=10 -9 ε=10 -8 ε=10 -7 ε=10 -6 ε=10 -5 ε=10 -4 ε=10 -3 ε=10 -2 ε=10 -1
Adaptive MR thresholding ε t ( We use the third order FV scheme designed in [START_REF] Aymard | A numerical method for transport equations with discontinuous flux functions: Application to mathematical modeling of cell dynamics[END_REF] and we couple it with the MR strategy on four levels. The domain consists of eight cell cycles with ten meshes in each cycle at the coarsest level in the x direction, and ten meshes in the y direction. The CFL number ensuring the stability of the scheme is set to 0.4. Note that the velocity coefficients depend on the solution (see equations (39-44) in Appendix A). The time step value is therefore time-dependent and computed at each time iteration. We adapt the stability rule designed in [START_REF] Aymard | Numerical simulation of the selection process of the ovarian follicles[END_REF] to the adaptive case so that

∆t n := min CF L min D ,(i,j) ∈ Sn ∆x max(|g n ,(i,j) |, |h n ,(i,j) |) , 1 max D ,(i,j) ∈ Sn |Λ n ,(i,j) | (36)
where g n ,(i,j) (respectively h n ,(i,j) ) is the discrete value of the speed g given by (39) (resp. h given by (40)) and Λ n ,(i,j)

is the discrete value of the source term Λ(x, y, U (t)) given by (41) on the mesh D ,(i,j) ∈ Sn . A simpler condition would consist in using the uniform case definition with the space step corresponding to the finest MR level

∆t n := min CF L ∆x L-1 max D ,(i,j) ∈ Sn (|g n ,(i,j) |, |h n ,(i,j) |) , 1 max D ,(i,j) ∈ Sn |Λ n ,(i,j) |
but since the speed and source coefficients are space dependent, the mesh has to be fully surveyed at each time iteration in any cases and formula (36) is more accurate and possibly less stringent. The numerical value of the biological constants appearing in definitions (39), ( 40) and ( 44) are gathered in Table 2. In the following paragraph we detail a set of numerical experiments which highlight the robustness of our numerical method in this biological set-up. We show in particular, that for this standard set of parameters spurious noise will appear at late times of the simulation unless the adaptive strategies presented earlier are used.

Numerical performances

The reference solution is computed on the finest level of resolution, comprising 51200 meshes and it takes 540 seconds on a mono processor Apple laptop to reach the final time t = 20. Snapshots of the solution at representative times are displayed in Figure 13, along with the adaptive mesh corresponding to ε = 0.01. Noticeably, at initial time (top panel), due to the periodicity of the MR reconstruction operator (see Eq. ( 22)), the mesh is refined near the right edge of the domain although the density is zero in this area.

As in the simple test case, we perform simulations with or without time adaptation of the MR threshold ( 19) and (34), and with or without the adaptive doubling strategy (20). Figure 14 dedicated to snapshots of the solution at t = 20 shows the influence of adaptive thresholds on the spurious noise removal. A zoom of the density computed with the different methods is displayed in the left panels. Both solutions look very similar although if we let the simulation run longer the solution in the bottom panel eventually deteriorates. In the right panels we display the adaptive grid superimposed to the density with a color code adapted to the amplitude of the spurious noise.

We then compute the L1 norm of the relative error (35) as a function of ε at different times, corresponding to the snapshots displayed in Figures 13 and14, with respect to the solution obtained using the uniform finest grid.

In Figure 15 we display the error curves against the MR threshold ε (left panel), the CPU gain (center panel) and the memory gain (right panel). The same legend as in Figure 9 is used, and described in Table 1. As in the simple 1D test case, the asymptotic behavior of the error in O(ε) is verified, except for extreme values of the threshold parameter ε.

In the left panels, the error curves corresponding to simulations with the adaptive doubling strategy (20) (with red +, blue * and cyan ) lie lower than their counterparts without the adaptive doubling strategy (with green ×, pink and black •). In the center and right panels, the performances of the different algorithms in terms of CPU time and memory requirements are displayed for a given accuracy measured by the relative error (35). The differences between the different strategies are not as spectacular as is the simple test case. Due to the differentiation phenomenon and the associated cell cycle exit, the overall proliferation rate is around 37 at the final time t = 20, much smaller -by a factor of 200-than in the 1D test case, so that the spurious noise is less important. Nonetheless, for late times, the curves corresponding to the simulation without MR nor doubling threshold adaptation (with green ×) clearly exhibit worse performances than the others. It is also noticeable that as time increases, the range of memory gain shifts to the left. Indeed, the size of the adaptive grid increases, because the density spreads out when it crosses the boundary between the proliferation and differentiation phase. The range of CPU gain also shifts accordingly to the left, since the computing time is directly and almost linearly related to the size of the grid.

The relation between the MR threshold ε and the adaptive doubling threshold δ is studied in Figure 16, where the relative error with the uniform grid solution is displayed against δ for different values of ε. Here again an optimal value is clearly visible. The curves in the right panel, where the product εE(δ, ε) is displayed against δ/ε are almost all superimposed, except for the two extreme values of ε. This corroborates the asymptotic behavior of the MR error in O(ε), which is always verified except if ε is too large (= 0.1) or too small (in that case the numerical accuracy of the computer penalizes the asymptotic behavior).

Comparison of the numerical performances for both models of mitosis

As previously explained in the Introduction and the beginning of Paragraph 4, the mitosis is an almost punctual event in time. Nevertheless cell proliferation is commonly modeled by a distributed gain term [START_REF] Perthame | Transport Equations In Biology[END_REF]. In this paragraph we study the influence of treating the mitosis in either a localized or distributed manner on the convergence of the MR algorithm with ε. In contrast with the simple model illustrated in Figure 7, we are now in a realistic biological context including the 2D dependence of the cell density on age and maturity variables. However, to make the comparison meaningful, the aging velocity is left constant in time and piecewise constant in age and maturity instead of using its general closed-loop model formulation (39). We can compute the equivalent distributed linear source term B, which results in doubling the cell mass after one cell cycle has elapsed

g(x, y, u) = γ 2 = 1.2 for (x, y) ∈ Θ 1 , 1 for (x, y) ∈ Θ 2 ∪ Ω 3 , (37)
B = 2 ln(2)γ 2 1 + γ 2 .
In order to reach a realistic mass gain before all cells have left the proliferation stage, we also reduce the time constant of the maturation velocity τ h = 0.2. Figure 17 displays the dependence of the relative error between the adaptive and uniform solutions with respect to the threshold ε, the CPU and the memory requirement. The left panels in Figure 17 show that the theoretical behavior in O(ε) is observed except for very small values of ε where the error estimates reaches the numerical accuracy. The reference solution on the uniform finest grid with 51200 meshes requires 535 seconds of CPU time. The center and right panels show that an average gain of 10 in CPU can be achieved, for a relative error of 10 -3 , while reducing the memory by a factor of six. Comparison of the top and bottom panels corresponding respectively to the localized and distributed models shows that the adaptive scheme exhibits comparable characteristics in both cases.

Conclusion

We have adapted a MR procedure to FV in order to simulate numerically non conservative transport equations with specific discontinuous flux conditions on interfaces between subdomains. In the absence of any specific care, we almost systematically observe the appearance of spurious numerical noise in long time simulations.

In this paper we have explained on a simple but generic model the reasons underlying this artifact and proposed a strategy to avoid it. We show that the doubling condition should be applied according to an adaptive rule only when the solution is above some threshold δ, to be taken of the order of the MR threshold ε. This stops the spurious noise appearance, while preserving all the robust characteristics of the MR method, namely its O(ε) asymptotic behavior, and interesting gains in CPU time and memory requirements. Accessorily, we have also qualitatively compared the modeling of biological mitosis by a discontinuous flux condition, to the distributed birth term often encountered in the literature, and shown that the later is also well handled by our FVMR scheme. This method will next be applied to real case simulations involving large dimension vector systems. 35) with respect to the threshold parameter ε (left panels), CPU gain (center panels) and memory gain (right panels). Top panels: Localized doubling with uncontrolled aging velocity (37), bottom panels: equivalent distributed doubling. The complete legend code is described in Table 1.

B Example of spurious numerical noise appearance

We show on a simple example that spurious noise can appear due to doubling transitions and requires refinement on the finest model, while the real discontinuities get smoothed out by the numerical scheme diffusion. Suppose for instance that at time t n the solution is locally almost piecewise constant, with a small discontinuity in the mesh D 0,K+1 requiring both this mesh and its immediate neighbors to be subdivided, as in Figure 18. To fix ideas, denote by α = 4ε/3 the size of the discontinuity

φ n 0,K-i = a, for i ≥ 1 φ n 1,2K = φ n 1,2K+1 = φ n 1,2K+2 = a φ n 1,2K+3 = φ n 1,2K+4 = φ n 1,2K+5 = a + α φ n 0,K+i = a + α, for i ≥ 3
The details at time n can be computed We have tested these details against the level dependent threshold ( 14) and we have found justifying the discretization on the coarsest and intermediate levels. Let us now update the solution at the next time step, using λ/2 = ∆t/∆x 1 on the intermediate level and λ/4 = ∆t/∆x 0 on the coarsest one (note that to update the solution on mesh (1, 2K) we will need to reconstruct the solution on its neighbor φn 1,2K-1 = a). We obtain the following values, also sketched out in the top panel of Figure 19 for instance say a = 4ε we have

|d n 1,K | = |d n 1,K+2 | < ε 1 , ε 1 ≤ |d n 1,K+1 | < ε 2 ,
                     φ n+1 0,K-i = φ n+1 0,K-1 = φ n 0,K-1 -λ 4 (φ n 0,K-1 -φ n 0,K-2 ) = a φ n+1 1,2K = φ n 1,2K -λ 2 (φ n 1,2K -2φ n 1,2K-1 ) = a + λa φ n+1 1,2K+1 = φ n 1,2K+1 -λ 2 (φ n 1,2K+1 -φ n 1,2K ) = a φ
   ε 2 ≤ |d n+1 1,K |, |d n+1 1,K+1 | < ε 1 , |d n+1 1,K+2 | < ε 1 , (49)
meaning that meshes D 1,2K and D 1,2K+1 have to be subdivided on the finest level while subdivisions D 0,K+1 and D 0,K+2 have to be aggregated again. The resulting reconstructed solution, obtained after applying Algorithms 1 and 2, is sketched out in bottom panel of Figure 19. This is of course a theoretical situation but it shows nevertheless how some spurious structure requiring fine meshes may appear in an upstream region and be well detached from the support of the genuine solution.

Figure 2 :

 2 Figure 2: Snapshots of the solution φ(t, x, y) at t = 20 on the whole computational domain and color code from 0 to 10 7 . Panels from top to bottom: a) Zoom on the support of the solution, MR with Haar reconstruction, b) MR with quadratic reconstruction, c) MR with quadratic reconstruction and doubling threshold strategy (20) d) finer grid no MR, e) full finer grid no MR. Notice the presence of the spurious noise (refined meshes upstream of and contiguous to the support of the solution in panel a) and b), as well as isolated upstream refined meshes in panel b)), and the absence of noise (in panel c)).

Figure 2 showsFigure 3 :

 23 Figure 3: Zoom on the support of the solution at t = 20 and computational grid with color code from 0 to to 10. From top to bottom a) MR with Haar reconstruction, b) MR with quadratic reconstruction, c) MR with quadratic reconstruction and doubling threshold (20). The thresholding error is more visible in the Haar case (panel a)) than in the quadratic reconstruction (panels b) and c)).

Figure 4 :

 4 Figure 4: Set-up for the p th interface. The characteristic curves in the right part of the domain transport either the initial condition if they cross the x-axis or the trace of the solution on the px s interface, multiplied by the k p coefficient, if they cross the vertical x = px s axis

Figure 5 Figure 5 :

 55 Figure 5: Mass with respect to time computed with the recursive definition (29), without and with the doubling strategy (20) for different values of the parameter δ. Full time range on the left and zoom on the final region on the right.

Figure 6

 6 Figure 6 displays the L1 norm over the time interval of the error between the thresholded and exact masses, normalized by the L1 norm of the mass of the exact solution

Figure 6 :

 6 Figure 6: L1 norm error (32) of mass with respect to the parameter δ of the doubling strategy (20) with (E(δ t )) and without (E(δ)) time adaptivity (33). The curve O(δ) indicates the theoretical asymptotic behavior which is a straight line of slope one in logarithmic axes.

Figure 7 :

 7 Figure 7: Mass over eleven cycles as a function of time for localized and distributed doubling. Third order scheme, four MR levels, ε = 10 -3 . The exact solution (29) is superimposed on the FVMR solution.

8 Figure 8 :

 88 Figure 8: Snapshot of the density φ(t, x, y) at time t = 7.48 (top panels) and t = 8 (bottom panels). Left panels: localized doubling, right panels : distributed gain over the cell cycle. The color bar giving the amplitude of the density (between 0 and 7000) is the same for all panels. The initial solution, given by (31), has a unit mass.

Figure 9 :

 9 Figure9: L1 norm of relative error (35) at different times, with respect to the MR threshold ε (left panel), the CPU gain (center panel) and the memory gain (right panel). Effect of adaptive MR threshold ε t and doubling threshold δ in the case of the simple 1D test case, 4 MR levels, third order FV scheme. The line labelled O(ε) in the left panel graphs indicates the theoretical asymptotic behavior of the error as a function of ε, which is a straight line of slope one in logarithmic axes. The complete legend code is described in Table1.

Figure 10 :

 10 Figure 10: Relation between the MR threshold parameter ε and the threshold δ in the adaptive doubling interface (20). L 1 norm of the error with solution computed on the uniform finest grid at final time. Simple 1D test case (1) with first order FV scheme. Left panels: error as a function of δ. Right panels: error/ε as a function of δ/ε, for different numbers of MR levels (top panels: L = 2, middle panels: L = 3, bottom panels: L = 4).

Figure 11 :Figure 12 :

 1112 Figure 11: Relation between the MR threshold parameter ε and the threshold δ in the adaptive doubling interface (20). L 1 norm of the error with the solution computed on the uniform finest grid at final time. Simple 1D test case (1), four MR levels, third order FV scheme. Left panels: error as a function of δ. Right panels: error/ε as a function of δ/ε. Top panels: constant MR threshold ε, bottom panels: adaptive MR threshold ε t .

Figure 13 :Figure 14 :

 1314 Figure 13: Snapshots of the density of granulosa cells computed with the FVMR strategy, adaptive doubling (20) and MR threshold[START_REF] Tomura | Contrasting quiescent G0 phase with mitotic cell cycling in the mouse immune system[END_REF], four MR levels, third order scheme, ε = 10 -2 , and biological constants from Table2. For each time, the adaptive grid is shown below the density, displayed with a time-dependent color code. The boundaries between biological cell cycles are marked with white lines.

Figure 15 :

 15 Figure 15: L1 norm of relative error (35) at different times, with respect to the MR threshold ε (left panel), the CPU gain (center panel) and the memory gain (right panel). Effect of adaptive MR threshold ε t and doubling threshold δ in the case of the biological model. The line labeled O(ε) in the left panels indicates the theoretical asymptotic behavior of the error as a function of ε which is a straight line of slope one in logarithmic axes. The complete legend code is described in Table1.

  Figure 15: L1 norm of relative error (35) at different times, with respect to the MR threshold ε (left panel), the CPU gain (center panel) and the memory gain (right panel). Effect of adaptive MR threshold ε t and doubling threshold δ in the case of the biological model. The line labeled O(ε) in the left panels indicates the theoretical asymptotic behavior of the error as a function of ε which is a straight line of slope one in logarithmic axes. The complete legend code is described in Table1.

1 Figure 16 :

 116 Figure 16: Relation between the MR threshold parameter ε and the threshold δ in the adaptive doubling interface (20). L 1 norm of the error with the solution computed on the uniform finest grid at final time. Biological model, four MR levels, third order FV scheme. Error as a function of δ (left), error as a function of δ/ε, error/ε as a function of δ/ε, for different values of ε. Constant MR threshold ε.

Figure 17 :

 17 Figure 17: L1 norm of relative error (35) with respect to the threshold parameter ε (left panels), CPU gain (center panels) and memory gain (right panels). Top panels: Localized doubling with uncontrolled aging velocity (37), bottom panels: equivalent distributed doubling. The complete legend code is described in Table1.

Figure 18 :

 18 Figure18: Set-up for Appendix B. Three level discretization with doubling interface located between meshes D 0,K-1 and D 0,K on the coarsest level. At time t n the solution is piecewise constant with a discontinuity of size α located between meshes D 1,2K+2 and D 1,2K+3 , resulting in a gradually adaptive grid.

Figure 19 :

 19 Figure 19: Discrete solution in example B. Top panel: updated solution at time t n+1 with formula (48). Bottom panel: solution computed with detail thresholding (49) according to Algorithms 1 and 2.

n+1 1 ,

 1 2K+2 = φ n 1,2K+2 -λ 2 (φ n 1,2K+2 -φ n 1,2K+1 ) = a φ n+1 1,2K+3 = φ n 1,2K+3 -λ 2 (φ n 1,2K+3 -φ n 1,2K+2 ) = a + α(1 -λ 2 ) φ n+1 1,2K+5 = φ n+1 1,2K+4 = φ n 1,2K+4 -λ 2 (φ n 1,2K+4 -φ n 1,2K+3 ) = a + α φ n+1 0,K+i = φ n+1 0,K+3 = φ n 0,K+3 -λ 4 (φ n 0,K+3 -φ n 0,K+2 ) = a + α(48)From these updated values we compute the averaged values on the coarsest level

Table 2 :

 2 Values of the parameters for the biological model simulations

	Parameter	Description	Value
		FSH plasma level (eq. (43))
	U min	minimum level	0.075
	U max	minimum level	0.15
	c	slope parameter	2.0
	M	abscissa of the inflection point	4.5
		Apoptosis source term (eq. (41))
	Λ	intensity factor	1
	Λ w	half width	0.05
	γ	scaling factor	0.1
	y s	cellular maturity threshold	0.3
		intrafollicular FSH level (eq. (44))
	b 1	basal level	0.08
	b 2	exponential rate	2.2482
	b 3	scaling factor	1450.
		Aging function (eq. (39))
	γ 1	rate	0.2
	γ 2	origin	1.0
	c 1		0.68
	c 2 ū	Maturation function (eq. (40))	0.08 0.02
	τ h		0.3

  If we use λ = 1/2 and the value α = 4ε/3 we obtain

	d n+1 1,K+2 = φ n+1 1,2K+4 -				1 8	φ n+1 0,K+1 + φ n+1 0,K+2 -	1 8	φ n+1 0,K+3
	= a + α -	1 8	(a +	α 2	λ 2	) + a + α -	1 8	(a + α) =	-αλ 32	+	α 8
				d n+1 1,K d n+1 1,K+1 = =	-a 32 a 4 +	ε 48
				d n+1 1,K+2 =	7ε 48
	we see that as soon as
				4	47 48	ε ≤ a < 16ε
	= a + λa -			1 8	a + a + a	λ 2	-	1 8	(a +	α 2	λ 2	) =	λa 2	+	αλ 32
	d n+1 1,K+1 = φ n+1 1,2K+2 -	1 8	φ n+1 0,K + φ n+1 0,K+1 -	1 8	φ n+1 0,K+2
	= a -	1 8	(a + a	λ 2	) + a +	α 2	λ 2	-	1 8	(a + α) = -	aλ 16	-	αλ 4	+	α 8

A Description of the model for follicular development

This paragraph presents the details of the biological model illustrated by the numerical simulations commented in Paragraph 4. It is meant to make the paper self-contained and introduce the notations. Details on the biological bases of the model can be found in the founding paper [START_REF] Echenim | Multi-scale modeling of the follicle selection process in the ovary[END_REF] and the review [START_REF] Clément | Multiscale modelling of ovarian follicular selection[END_REF]. Since our goal here is to design the coupling of MR with adaptive doubling strategy, the vector nature of the model is not relevant, so that we only present its scalar version, referring the reader to [START_REF] Aymard | Numerical simulation of the selection process of the ovarian follicles[END_REF][START_REF] Aymard | A numerical method for transport equations with discontinuous flux functions: Application to mathematical modeling of cell dynamics[END_REF] for a more general presentation of the numerical scheme. Let us denote φ(t, x, y) the density of granulosa cells at time t, age x and cell maturity y. It satisfies the following equation :

set in the computing domain Ω in the (x, y) plane,

where N c is the number of cell cycles and D c is the duration of one cycle. The domain Ω is divided in 2N c + 1 zones:

and Ω 3 , corresponding to different cell states illustrated in Figure 12 and hence different definitions of the speeds and source terms. Phase Ω p 2 in the model aggregates the three latest phases (S, G2, M) of the

Each cell cycle consists of the Ω p 1 ∪ Ω p 2 subdomain and Θ i for i = 1, 2 denotes the disconnected union of the N c corresponding phases Ω i p , for p = 1, . . . , N c . The aging function g appearing in (38) is defined by

where γ 1 , γ 2 are real positive constants. The maturation function h is defined by

where τ h , c 1 , c 2 and ū are real positive constants.

The source term, that represents cell loss through apoptosis, is defined by

where Λ, y s and γ are real positive constants, and χ E (y) is the characteristic function of subset E.

The equations in the PDE system (38) are linked together through the argument u(t) appearing in the speeds g(x, y, u) and h(x, y, u) and the argument U (t) appearing in the source term Λ(x, y, U ). U (t) and u(t) represent respectively the plasma FSH level and the locally bioavailable FSH level and depend on the first maturity moment of the density

The plasma FSH level U (t) showing up in the arguments of the source term in (38) is defined by

where U min , c and M are real positive constants. The locally bioavailable FSH level u(t) showing up in the arguments of the speeds in (38) is defined by

where b 1 , b 2 and b 3 are real positive constants. The precise definition of the required transmission conditions along the successive cell cycles of the domain has been formulated in [START_REF] Shang | Cauchy problem for multiscale conservation laws : Applications to structured cell populations[END_REF]. For each cycle p = 1, . . . , N c ,

• the flux on the x-axis is continuous on the interface between Ω p 1 and

• The flux is doubling on the interface between Ω p 2 and Ω p+1 1 , which accounts for the birth of two daughter cells from one mother cell at the end of each cell cycle (γ 1 u + γ 2 )φ(t, x + , y) = 2φ(t, x -, y), x = pD c , 0 ≤ y ≤ y s .

(46)

• A homogeneous Dirichlet condition holds to the north of the interface between Ω p 2 and Ω 3 φ(t, x, y + s ) = 0, (p -1/2)D c ≤ x ≤ pD c .

(47)