Benjamin Aymard

Frédérique Clément

Frédéric Coquel

Marie Postel

Aymard Benjamin

Clément Frédérique

Coquel Frédéric

Postel Marie

Numerical simulation of the selection process of the ovarian follicles

This paper presents the design and implementation of a numerical method to simulate a multiscale model describing the selection process in ovarian follicles. The PDE model consists in a quasi-linear hyperbolic system of large size, namely N f × N f , ruling the time evolution of the cell density functions of N f follicles (in practice N f is of the order of a few to twenty). These equations are weakly coupled through the sum of the first order moments of the density functions. The timedependent equations make use of two structuring variables, age and maturity, which play the roles of space variables. The problem is naturally set over a compact domain of R 2 . The formulation of the time-dependent controlled transport coefficients accounts for available biological knowledge on follicular cell kinetics. We introduce a dedicated numerical scheme that is amenable to parallelization, by taking advantage of the weak coupling. Numerical illustrations assess the relevance of the proposed method both in term of accuracy and HPC achievements.

Résumé. Ce document présente la conception et l'implémentation d'une méthode numérique servant à simuler un modèle multiéchelle décrivant le processus de sélection des follicules ovariens. Le modèle EDP consiste en un système hyperbolique quasi linéaire de grande taille, typiquement N f × N f , gouvernant l'évolution des fonctions de densité cellulaire pour N f follicules (en pratique N f est de l'ordre de quelque uns à une vingtaine). Ces équations d'évolution utilisent deux variables structurantes, l'âge et la maturité, qui jouent le rôle de variables d'espace. Le problème est naturellement posé sur un domaine compact de R 2 . La formulation du transport à coefficients variables au cours du temps en fonction du contrôle est issue des connaissances disponibles sur la cinétique cellulaire au sein des follicules ovariens. Nous présentons un schéma numérique dédié au problème parallélisable en tirant parti du faible couplage. Des simulations numériques permettent d'évaluer la pertinence de la méthode proposée tant en terme de précision que de calcul haute performance.

Introduction

This work is motivated by a mathematical modelling approach of a complex physiological system, the development of ovarian follicles. The model describes both the cell dynamics within each follicle and the competition process within the population of follicles. The resulting model ( [START_REF] Echenim | Multi-scale modeling of the follicle selection process in the ovary[END_REF]) is a large scale system of weakly coupled quasi-linear transport equations, where integro-differential terms occur both in the velocity and source term. The coupling terms account for the endocrine-based dependence of one follicle dynamics on all other developing 1. Biological and biomathematical background

Mesoscopic scale: ovarian follicles

Ovarian follicles are spheroidal tissular structures sheltering the oocytes (figure 1). Each follicle goes through a development process, composed of two main parts : the basal development and the terminal development (figure 2). During basal development, the follicles are independent of hormonal support from the pituitary gland, while, during the terminal part of their development, they become dependent on FSH (follicle stimulating hormone) supply. Follicle development can end up either by ovulation, in the best case, or much more often by degeneration through atresia. The commitment of a follicle to either ovulation or atresia is driven by the changes occurring in the follicular cell population. After exiting the pool of quiescent primordial follicles, ovarian follicles enter a several-month long process of growth and maturation, that either ends up by ovulation or degeneration through atresia. In the basal part of development, follicle growth is mainly due to the oocyte enlargement and increase in the number of granulosa cell layers. In the final part of development, granulosa cells progressively stop proliferating and follicle growth is mainly due to the enlargement of the antrum, a liquidfilled cavity in the center of the follicle.

Macroscopic scale : competition process and endocrine feedback loop

In some sense, the terminal part of follicle development can be considered as a competition between follicles for FSH resource. FSH controls follicle development, while FSH levels are in turn controlled by two hormones secreted by the ovary, inhibin and estradiol, whose plasma levels are determined by the summed contributions of all maturing follicles (figure 1). This ovarian hormonal feedback induces a drop in FSH levels that will penalize the follicles, except those (in the poly-ovulating situations) or that (in the mono-ovulating ones) that are sufficiently mature to survive in a FSH poor environment. The rising levels of estradiol finally triggers the ovulatory surge that leads to the ovulation of the surviving follicles.

Microscopic scale : granulosa cell kinetics

Each granulosa cell can be encountered in either of three different cell states : proliferation, differentiation or apoptosis (programmed cell death) (figure 3). At the beginning of terminal development, most granulosa cells are progressing along the cell division cycle, that can be split into a G1 phase, where cells are sensitive to FSH control, and a SM phase, where cells are preparing for mitosis and are insensitive to FSH control. At the end of mitosis, one single mother cell gives birth to two daughter cells. During the FSH sensitive phase, cells become more and more mature, up to a point where they reach a threshold maturity and exit the cell division cycle. At this time they are exposed to a great risk of apoptosis, and may die if the FSH environment is not favorable. After exiting the cell cycle, the cells stop proliferating definitively, but their maturity still increases, so that they contribute more and more to hormone (and especially estradiol) secretion.

Biomathematical model

Φ(a, γ, t) = (φ 1 (a, γ, t), ..., φ N f (a, γ, t)),
and the computing domain Ω in the (a, γ) plane,

Ω = {(a, γ), 0 ≤ a ≤ N c × D a , 0 ≤ γ ≤ 1}
where N c is the number of cell cycles and D a is the duration of one cycle. The different cell states described in paragraph 1.1.3 take place in three subdomains as illustrated in Figure 3

     G1 = {(a, γ) ∈ Ω, pD a ≤ a ≤ (p + 1/2)D a , p = 0, . . . , N c -1, 0 ≤ γ ≤ γ s }, SM = {(a, γ) ∈ Ω, (p + 1/2)D a ≤ a ≤ (p + 1)D a , p = 0, . . . , N c -1, 0 ≤ γ ≤ γ s }, D = {(a, γ) ∈ Ω, γ s ≤ γ}.

Hyperbolic system

If we consider, for each follicle, a conservation law for the granulosa cell population density, we obtain the following hyperbolic system, satisfied by Φ:

                           ∂φ 1 (a, γ, t) ∂t + ∂(g(a, γ, u 1 (t))φ 1 (a, γ, t)) ∂a + ∂(h(a, γ, u 1 (t))φ 1 (a, γ, t)) ∂ γ = -λ(a, γ, U (t))φ 1 (a, γ, t) . . . ∂φ f (a, γ, t) ∂t + ∂(g(a, γ, u f (t))φ f (a, γ, t)) ∂a + ∂(h(a, γ, u f (t))φ f (a, γ, t)) ∂ γ = -λ(a, γ, U (t))φ f (a, γ, t) . . . ∂φ N f (a, γ, t) ∂t + ∂(g(a, γ, u N f (t))φ N f (a, γ, t)) ∂a + ∂(h(a, γ, u N f (t))φ N f (a, γ, t)) ∂ γ = -λ(a, γ, U (t))φ N f (a, γ, t) (1)
The initial condition is given by the distribution of the granulosa cell populations at the initial time

Φ(a, γ, 0) = Φ 0 (a, γ)
where Φ 0 is compactly supported in ]0, N c D a [×]0, 1[. Since the functional domain Ω is compact, we need to express boundary conditions. Considering some qualitative properties of the maturation function h, defined below, the horizontal top boundary γ = 1 is never reached. Similarly, we can fix the number of cell cycles N c in accordance with the aging function g, so that the right vertical boundary a = N c D a is never reached. For sake of computing simplicity, we therefore assume spatial periodicity on the outer boundaries of Ω in the numerical simulations. The aging function g appearing in (1) is defined by

g(a, γ, u) = g 1 u + g 2 for (a, γ) ∈ G1 1 for (a, γ) ∈ SM ∪ D (2)
where g 1 , g 2 are real positive constants. The maturation function h is defined by

h(a, γ, u) = (-γ 2 + (c 1 γ + c 2 )(1 -exp( -u ū ))) for (a, γ) ∈ G1 ∪ D 0 for (a, γ) ∈ SM (3)
where c 1 , c 2 and ū are real positive constants. The source term, that represents cell loss through apoptosis, is defined by

λ(a, γ, U ) =    K exp(-( (γ -γ s ) 2 γ )) × (1 -U ) for (a, γ) ∈ G1 ∪ D 0 for (a, γ) ∈ SM (4)
where K, γ s and γ are real positive constants.

Remark. The precise values of the constants depends on species and breeds; they will be fixed later.

Closure equations

The equations in the PDE system (1) are linked together through the arguments u f (t) and U (t) appearing in the speeds g(a, γ, u) and h(a, γ, u) and in the source term λ(a, γ, u). U (t) and u f (t) represent respectively the plasma FSH level and the locally bioavailable FSH level and depend on some maturation moments of the densities.

• follicular cell mass (local, one by follicle)

m 0 (f, t) = 1 0 NcDa 0 φ f (a, γ, t)dadγ, 1 ≤ f ≤ N f ( 5
)
This moment is an interesting quantity since it corresponds to an observable variable (the total cell number in one follicle) even if it does not come as such into play in the equations. • follicular maturity (local, one by follicle)

m(f, t) = 1 0 NcDa 0 γφ f (a, γ, t)dadγ, 1 ≤ f ≤ N f (6)
• ovarian maturity (global, shared by all follicles)

M (t) = N f f =1 m(f, t), (7)
• plasma FSH level (global, shared by all follicles)

U (t) = U min + 1 -U min 1 + exp(c(M (t) -m)) , (8)
where U min , c and m are real positive constants

• locally bioavailable FSH level u f (t) = min b 1 + e b2m(f,t) b 3 , 1 U (t), (9)
where b 1 , b 2 and b 3 are real positive constants.

It is important to notice that the aging function ( 2) is discontinuous (in general) on the interfaces G1 -SM , and that the maturation function ( 3) is discontinuous (in general) on the interfaces SM -D (Figure 3). It is necessary to introduce transmission conditions in order to overcome possible failure of uniqueness due to the discontinuities in coefficients (see [START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. I. The scalar case[END_REF] or [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF] for instance). The precise definition of the required transmission conditions has been addressed in the paper by Peipei Shang ( [17]). We suppose that, for each cycle p = 1, . . . , N c and each follicle f = 1, . . . , N f , the flux on the a-axis is continuous between the phases G1 and SM

φ f (t, a + , γ) = (g 1 u f + g 2 )φ f (t, a -, γ), a = (p -1/2)D a , 0 ≤ γ ≤ γ s , (10)
and that the flux is doubling on the interfaces SM -G1, which accounts for the birth of new cells at the end of each cell cycle

(g 1 u f + g 2 )φ f (t, a + , γ) = 2φ f (t, a -, γ), a = pD a , 0 ≤ γ ≤ γ s . (11)
Finally, we suppose homogeneous Dirichlet condition to the north of the interface

SM -D φ f (t, a, γ + s ) = 0, (p -1/2)D a ≤ a ≤ pD a . (12)
The time dependent ovarian maturity M (t) defined in ( 7) is compared to the ovarian maturity threshold, denoted by M o , to define the final time, beyond which the competition is over

t f inal = inf{t, M (t) ≥ M o }. (13)
The follicles are then sorted into two classes : ovulatory if the maturity is higher than the follicle maturity threshold, denoted by M f , atretic otherwise.

Numerical Method

Discretization

In this paragraph the cycle duration is set to D a = 1. We denote by ∆a (respectively ∆γ ) the space step in the age (respectively maturity) direction. In practice we choose ∆a = ∆γ. The discretization step ∆γ and the cellular maturity threshold γ s must be chosen so that the interfaces where the speed coefficients (2,3) are discontinuous fall on grid points covering the computational domain of size [N c , 1]. Denoting by N m the number of grid cells by half granulosa cell cycle, we set N γ = 2N m and ∆a = ∆γ = 1/N γ and we introduce the dedicated2 notations for grid points (a k , γ l ) and mesh centers (a k+1/2 , γ l+1/2 )

a k = k∆a, a k+1/2 = (k + 1/2)∆a, for k = 0, . . . , N c × N γ , (14)
γ l = l∆γ, γ l+1/2 = (l + 1/2)∆γ, for l = 0, . . . , N γ . ( 15
)
Considering that the time step ∆t n may change at every iteration, in order to preserve stability, the time discretization is defined by

t 0 = 0, t n+1 = t n + ∆t n , for n = 0, . . . , N t (16)
with N t such that t Nt = t final . The unknowns are the approximate mean values of the density vector in each grid mesh

Φ n k,l ≈ 1 ∆a∆γ a k+1 a k γ l+1 γ l
Φ(a, γ, t n )dγda, for k = 0, ..., N c N γ -1, and l = 0, ..., N γ -1 whose components φ n f,k,l are the discrete density values for each follicle.

Macroscopic scale : piecewise constant approximation of the hormonal control

We define the approximation of the control terms (5 -9) at each time step n = 0, . . . , N t

m n 0,f = ∆a∆γ Nγ -1 l=0 NcNγ -1 k=0 φ n f,k,l , for f = 1, . . . , N f , (17)
m n f = ∆a∆γ Nγ -1 l=0 γ l+1/2 NcNγ -1 k=0 φ n f,k,l , for f = 1, . . . , N f , (18)
M n = N f f =1 m n f , (19)
U n = U min + 1 -U min 1 + exp(c(M n -m)) , (20)
u n f = min(b 1 + e b2m n f b 3 , 1)U n , for f = 1, . . . , N f . (21)

Mesoscopic scale : finite volume scheme

We use a splitting strategy in order to compute the solution of the PDE system [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF], which amounts to a convective equation combined with a source equation, for each follicle f = 1, . . . , N f

∂ t φ f (a, γ, t) + ∂ a (g(a, γ, u f (t))φ f (a, γ, t)) + ∂ γ (h(a, γ, u f (t))φ f (a, γ, t)) = 0 (convective part), ( 22
)
∂ t φ f (a, γ, t) = -λ(a, γ, U (t))φ f (a, γ, t) (source part). ( 23
) 2.3.1. Convective part
The convection part ( 22) is treated with a classical finite volume method. The approximate mean values of the solution at t = 0 are initialized using a midpoint formula, accurate at the order 2 in space

Φ 0 k,l = Φ 0 (a k+1/2 , γ l+1/2 ).
Using the integral form of the conservation law

t n+1 t n a k+1 a k γ l+1 γ l (∂ t φ f (a, γ, t) + ∂ a (g(a, γ, u f (t))φ f (a, γ, t)) + ∂ γ (h(a, γ, u f (t))φ f (a, γ, t)))dγdadt = 0,
we obtain a recursion on the approximate density of each follicle, where we drop the index f for clarity sake

φ n+1 k,l = φ n k,l - ∆t ∆a G k+1,l+ 1 2 (Φ n ) -G k,l+ 1 2 (Φ n ) - ∆t ∆γ H k+ 1 2 ,l+1 (Φ n ) -H k+ 1 2 ,l (Φ n ) (24)
where G k,l+ 1 2 (respectively H k+ 1 2 ,l ) is the numerical flux across the vertical edge [(a k , γ l ), (a k , γ l+1 )] (respectively the horizontal edge [(a k , γ l ), (a k+1 , γ l )]) defined by

         G k,l+ 1 2 (Φ n ) ≈ 1 ∆t n ∆γ t n+1 t n γ l+1 γ l g(a k , γ, u f (t))φ(a k , γ, t)dγdt, H k+ 1 2 ,l (Φ n ) ≈ 1 ∆t n ∆a t n+1 t n a k+1 a k h(a, γ l , u f (t))φ(a, γ l , t)dadt. ( 25
)
The dependence of the flux functions on all the follicle densities through the control term u f is emphasized by the Φ n argument. It appears in the explicit in time approximations of the speeds (2), and (3) at the center of the mesh

g n k,l = g(a k+ 1 2 , γ l+ 1 2 , u n f ), h n k,l = h(a k+ 1 2 , γ l+ 1 2 , u n f ).
Note that the transmission conditions [START_REF] Frey | Medit, an interactive mesh visualization software[END_REF] and [START_REF] Godlewski | The numerical interface coupling of nonlinear hyperbolic systems of conservation laws. I. The scalar case[END_REF] on the interfaces where these speed coefficients (2) and ( 3) are discontinuous are exactly treated with the method described in the article of Godlewski and Raviart ( [11]). The numerical fluxes (25) are designed using a limiter strategy. Indeed, it is well known that first order schemes, like the Godunov scheme, are diffusive, and that second order schemes, like Lax Wendroff scheme, generate oscillations in the neighborhood of discontinuities. In order to get a stable as well as precise scheme, we take a weighting of a low order scheme and a high order scheme, and we define limited numerical fluxes

G = G Low + ℓ(r G )(G High -G Low ), H = H Low + ℓ(r H )(H High -H Low ), ( 26
)
where ℓ is a limiter function, for example van Leer function

ℓ(r) = r + |r| 1 + r and r G , r H are r G =              g k-1,l φ k-1,l -g k-2,l φ k-2,l g k,l φ k,l -g k-1,l φ k-1,l if g k-2
,l ≥ 0 and g k-1,l ≥ 0 and g k,l ≥ 0,

g k+1,l φ k+1,l -g k,l φ k,l g k,l φ k,l -g k-1,l φ k-1,l if g k-1
,l ≤ 0, and g k,l ≤ 0 and g k+1,l ≤ 0, 0 otherwise,

r H =              h k,l-1 φ k,l-1 -h k,l-2 φ k,l-2 h k,l φ k,l -h k,l-1 φ k,l-1 if h k,l-2 ≥ 0 and h k,l-1 ≥ 0 and h k,l ≥ 0, h k,l+1 φ k,l+1 -h k,l φ k,l h k,l φ k,l -h k,l-1 φ k,l-1 if h k,l-1 ≤ 0 and h k,l ≤ 0 and h k,l+1 ≤ 0, 0 otherwise.
These ratios are good indicators of the regularity of the function in each direction (see [START_REF] Sweby | High resolution schemes using flux limiters for hyperbolic conservation laws[END_REF]). In fact, a steep gradient or a discontinuity gives a ratio far from 1, whereas a smooth function gives a ratio close to 1.

The first order fluxes entering equation ( 26) are the Godunov fluxes

     G Low k,l+ 1 2 (Φ n ) = (g n k-1,l ) + φ n k-1,l + (g n k,l ) -φ n k,l , H Low k+ 1 2 ,l (Φ n ) = (h n k,l-1 ) + φ n k,l-1 + (h n k,l ) -φ n k,l ,
and the high order fluxes are the Lax Wendroff ones

       G High k,l+ 1 2 (Φ n ) = 1 2 (g n k-1,l φ n k-1,l + g n k,l φ n k,l ), H High k+ 1 2 ,l (Φ n ) = 1 2 (h n k,l-1 φ n k,l-1 + h n k,l φ n k,l ).
Consequently, since ∆a = ∆γ, the CFL condition which guarantees the stability is

∆t n f,F V S ≤ CF L ∆γ max k,l (|g k,l |, |h k,l |) (27)
with CF L ≤ 1 2 .

Source part

The source part (23) of the PDE system is explicitly dealt with

φ n+1 k,l = φ n k,l -∆tλ(a k+ 1 2 , γ l+ 1 2 , U n )φ n k,l , (28)
which implies a stability condition on the time step, which we strengthen to enforce positivity

∆t n f,source ≤ 1 max k,l |λ(a k+ 1 2 , γ l+ 1 2 , U n )| . ( 29
)

Macro/Meso scale Parallelization

In order to speed up the numerical simulations, the code is implemented on a parallel architecture. Every follicle follows its own dynamic independently of the others, except that it shares with them the FSH resource, and that it contributes to the tuning of plasma FSH level. Hence it is natural to use a SIMD (single instruction, multiple data) strategy, with one follicle by process. The only communications that have to be done at each time step concern the ovarian maturity, through a reduction operation (sum) and the current time step ∆ n , through a reduction operation (min). Each process computes the maximum time step satisfying both conditions (27) and ( 29)

∆t n f = min{∆t n f,F V S , ∆t n f,source } (30)
relevant for its own follicle, and the time step used by all processes is Denoting by E the evolution operator, that consists in steps ( 2) and (3), the second order in time is achieved by a second order Runge-Kutta method (Heun)

∆t n = min f {∆t n f }. (31
φ * = E(φ n ), φ * * = E(φ * ), φ n+1 = 1 2 (φ n + φ * * ).

Numerical simulations

The code was developed in C++ using the parallel library MPI. The tests were performed on the Jacques-Louis Lions Laboratory super calculator SGI Altix UV 100. The current configuration is 64-core 2 GHz. The visualizations were made with MEDIT [START_REF] Frey | Medit, an interactive mesh visualization software[END_REF] and gnuplot.

Initial condition

For the initial condition, we take a gaussian function centred in C = (C a , C γ )

φ 0 (a, γ) = 1 2πσ exp 1 2 (a -C a ) 2 σ 2 + (γ -C γ ) 2 σ 2 . ( 32
)
The variance σ 2 can be chosen so that φ 0 is smooth enough to be used for testing the convergence rate.

Sets of default parameters

Since one goal of this work was to improve the computational method, we used the same set of parameters as in [START_REF] Hombourger | Rapport de stage de 3 me année. Modélisation multi-échelle du développement folliculaire[END_REF] to be able to compare at least qualitatively with our previous simulations outputs. We distinguish two sets of parameters, the global ones, which are the same for all follicles (space discretization, CFL condition,. . . ) in Table 1 and the local ones (initial condition, velocity parameters,. . . ), which can depend on the follicle, in

.3. Time evolution of the cell density for one follicle

The behavior of the cell density of one single follicle as a function of time is first studied, in order to check some qualitative properties. The full movie can be downloaded from http://www.ljll.math.upmc.fr/aymard/film.gif and four snapshots are displayed on Figure 4. Note that on the snapshots the color code is time dependent while in the movie it is set once and for all at initial time. The spatial discretization for this test uses N m = 30 cells per half granulosa cell cycle, and the domain Ω contains five cycles. The grid size is therefore 5 × (2 × 30) 2 = 18000.

The final time defined by ( 13) is not reached and the simulation is stopped at t max = 4 after 570 time steps. This corresponds to the time required to cover four cell cycles. On the first snapshot 4 a), we can observe the initial condition, at the center C = (0.25, 0.15) of the first G1 phase. On the second snapshot 4 b), we see that the density has moved to the second cycle and that it has been split into a fully differentiated part and a still proliferating part, that lies in the SM phase. It is worth noticing that, consistently with the model, there is no crossing of the SM -D interface. On the third snapshot 4 c), we can see that the cell density is doubling at the end of the cell cycle, on the SM -G1 interface. On the last snapshot 4 d), at the end of the simulation, the density is concentrated in the D phase above the fifth cell cycle. Even if all cells have exited the cell cycle, we can distinguish three different density clouds, each of which being issued from one of the previous cycles (respectively the second third and fourth cycle from top to down). This behavior of the numerical solution is in accordance with the previous simulations presented in [START_REF] Hombourger | Rapport de stage de 3 me année. Modélisation multi-échelle du développement folliculaire[END_REF].

Competition between ten follicles

The second test consists in a more realistic simulation of a competition process between ten follicles. The parameters defining the plasma FSH level are set to U min = 0.9 and c = 10 and the intensity of the source term is set to K = 1. The follicles are distinguished by the parameter defining their aging rate (2) at the origin. We consider a range of g 2 values running from 0.5 to 0.95, with a 0.05 increment from one follicle to the other g 2 = 0.5, 0.55, . . . , 0.95.

The maximum time is set to t max = 1.5. For all other parameters we use the values in Tables 1 and2. Figure 5 represents a) the plasma FSH level U (t) defined by (8) b) the locally bioavailable FSH level u f (t) defined by (9) c) the cell mass m 0 (f, t) of each follicle defined by (5), d) the ovarian maturity M (t) defined by ( 7) and e) the maturity m(f, t) of each follicle defined by [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF]. Figure 5.d) shows that the ovarian maturity reaches the threshold defined by [START_REF] Hovhannisyan | On the stability of fully adaptive multiscale schemes for conservation laws using approximate flux and source reconstruction strategies[END_REF] with M o = 5 around t = 1.2. Figure 5.e) shows that at the final time, only three follicles have reached the ovulatory stage, where the follicle maturity defined by ( 6) is higher than M f = 0.7. They correspond to the values of the parameter g 2 = 0.5, 0.55 and 0.63 .

Convergence rate test

We now turn to the validation of the code, which consists in numerically verifying the asymptotic order of convergence when the time step ∆t and space discretization step ∆γ go to zero. We use the parameters in Tables 1 and2, except for the gaussian function, whose variance is set to σ 2 = 0.002, so that it is very smooth. The simulation is stopped at t max = 0.05, which allows us to reduce the number of cell cycles to N c = 1 and to discretize the solution on a square grid N γ × N γ . We compute the solution for six different levels discretizations We first study the convergence of the numerical scheme when the initial condition is centred in the SM phase C = (0.7, 0.15) and the final time is small enough for the density to remain in this phase. The maturity function (3) and the source term (4) are both null, and the aging function ( 2) is constant and equal to one. For such a constant linear transport, we can compute the exact solution at the final time t Nt = t max using the characteristic method as well as the error in L 1 -norm of the numerical solution

E(∆γ) = ∆γ 2 Nγ -1 k=0 Nγ -1 l=0 φ Nt k,l -φ 0 (a k+ 1 2 -t Nt , γ l+ 1 2 ) .
The error curve in log scale is superposed with the theoretical order O(∆γ 2 ) in Figure 6. ). Blue curve : L1 norm error between the exact solution and the solution with space step ∆γ.

Full model convergence test

We now set the initial condition near the center C = (0.2, 0.15) of the G1 phase, still with a maximum time t max = 0.05 ensuring that the cell density remains in this zone. In that case the source term (4) and the age and maturity speeds are no more constant, and we no longer know the exact solution of the PDE. The error is computed using the solution with N γ = 5120 as reference solution, that is compared to the solution on the other meshes. Since we have used a dyadic refinement, the size N ref of the reference mesh in one direction is The asymptotic order O(∆γ 1.95 ) which best fits the behavior of E(∆γ) is displayed in Figure 7. ) E

Figure 7. CV rate on the G1 phase. Green curve : ∆γ → O(∆γ 1.95 ). Blue curve : L1 norm error between the reference solution and the solution with a space step ∆γ.

HPC test

The improvement in terms of computing time provided by the parallelization is tested through a set of simulations involving an increasing number of identical follicles. The grid for these simulations consists in one cycle of 200 × 200 cells, and the simulation goes on for 900 times steps. In this experiment, we disposed of enough processors to do the computing with one follicle by processor and since the time command is used to monitor the computing time, the output of the program has been commented out. Figure 8 a) shows the real computing time, which is basically the time elapsed from the beginning of the computation. Figure 8 b) shows the user computing time, which cumulates all the processors computing time. The computing time cannot be shorter than that required for a computation involving only one, this constitutes the so-called theoretical limit. The fact that the real computing time remains close to the computing time for one follicle and that the user computing time grows linearly with the number of follicles is very encouraging. The communications between processors appear not to affect significantly the gain in computing time due to the parallel computing. A small increase in real computing time can be noticed beyond eight follicles. This is due to the super calculator architecture, where the processors are pooled in eight processor nodes. The communications within one node are faster than across different nodes. This leads to a threshold effect observed as soon as more than eight processors are needed, therefore involving more than one single node.

Conclusion

This paper summarizes preliminary works in the development of a dedicated software to illustrate numerically the development of follicles. The uniform grid numerical method have been successfully tested in terms of robustness, accuracy and scalability on parallel architecture. As part of a challenging project involving biologists mathematicians and computer scientists, the different situations engendered by the model (mono-ovulation, poly-ovulation or anovulation) from given combinations of parameters will now be systematically and intensively tested. In order to achieve this goal within realistic delay the HPC aspect of the method must be enriched using adaptive mesh refinement. The multiresolution method developed in [START_REF] Cohen | Fully adaptive multiresolution finite volume schemes for conservation laws[END_REF] and further tested and extended in [START_REF] Coquel | Entropy-satisfying relaxation method with large time-steps for Euler IBVPs[END_REF] or [START_REF] Hovhannisyan | On the stability of fully adaptive multiscale schemes for conservation laws using approximate flux and source reconstruction strategies[END_REF] is currently implemented in this new configuration, where the presence of source terms and integral terms will require special considerations as exemplified for instance in [START_REF] Chiavassa | Cost-effective multiresolutions schemes for shock computations[END_REF].

Figure 1 .

 1 Figure 1. Follicle development as a multiscale process. Left : endocrine feedback loop between the ovaries and the hypothalamo-pituitary-axis. Middle : schematic 2D view of an ovarian follicle. Right : different cell states encountered by the granulosa cells during follicular development.

Figure 2 .

 2 Figure2. Follicular development : from basal to terminal development. After exiting the pool of quiescent primordial follicles, ovarian follicles enter a several-month long process of growth and maturation, that either ends up by ovulation or degeneration through atresia. In the basal part of development, follicle growth is mainly due to the oocyte enlargement and increase in the number of granulosa cell layers. In the final part of development, granulosa cells progressively stop proliferating and follicle growth is mainly due to the enlargement of the antrum, a liquidfilled cavity in the center of the follicle.

Figure 3 .

 3 Figure 3. Division of the spatial domain according to the cell phases for unit cycle duration, D a = 1. Left : schematic view of a granulosa cell cycle. Right : (a, γ) plane. The bottom part represents the successive cell cycles, each composed of the G1 and SM phases. The top part corresponds to the differentiation phase.

3

 3

NFigure 4 .

 4 Figure 4. Snapshots of the cell density at different times, starting from a smooth initial condition in the G1 phase. The color code is time dependent. The full movie is available at http://www.ljll.math.upmc.fr/aymard/film.gif

Figure 5 .

 5 Figure 5. Competition between ten follicles differing by their aging function parameter g 2

Figure 6 .

 6 Figure 6. Convergence rate in the SM phase. Green curve : ∆γ → O(∆γ 2 ). Blue curve : L1 norm error between the exact solution and the solution with space step ∆γ.

  always a power of two times the size N ∆γ of the current mesh. Denoting by P = N ref N ∆γ the ratio between the current discretization and the reference finest one, we estimate the discretization error by E(∆γ) = ∆γ 2

Figure 8 .

 8 Figure 8. Parallelization performance test. Simulation of an increasing number of identical follicles. One follicle per processor. Green curve : theoretical limit. Blue curve : computing time (in seconds). Left : elapsed time from the start of the computation. Right : cumulative computed time on all processors

  ) 2.5. Parallel algorithm and improvement to order 2 in time

	Finally, we get the following parallel algorithm for each time step:
	(1) Compute time step
	• Compute local time step ∆t n f (30) • Communication : Compute the minimum of local time steps ∆t n (31)
	(2) Update solution
	• Flux computation (26),
	• Convection computation (24),
	• Source computation (28),
	(3) Update control terms
	• Compute local maturity (18),
	• Communication : compute the sum of local maturities (19),
	• Compute FSH plasma level (20),
	• Compute locally bioavailable FSH rate (21).

Table 2 .

 2

	Parameter	Description	Value
	N m	number of grid cells by half granulosa cell cycle
	N c	number of cycles	5
	CFL	CFL condition	0.45
	M f	follicular maturity threshold	3.0
	M o	ovarian maturity threshold	7.0
		FSH plasma level (eq. (8))
	U min	minimum level	0.075
	c	slope parameter	2.0
	m	abscissa of the inflection point	4.5
		Apoptosis source term (eq. (4))
	K	intensity factor	6.0
	γ	scaling factor	0.2
	t max	maximum time (to avoid excessively long computations)

Table 1 .

 1 Values of the global parameters used in Figures 4 to 8The number of follicles, N f , and the number of grid cells by half granulosa cell cycle, N m , may depend on the simulation.

	Parameter	Description	Value
	γ s	cellular maturity threshold	0.3
		intrafollicular FSH level (eq. (9))
	b 1	basal level	0.08
	b 2	exponential rate	2.25
	b 3	scaling factor	1450.
		Aging function (eq. (2))
	g 1	rate	4.0
	g 2	origin	0.7
		Maturation function (eq. (3))
	c 1		0.68
	c 2		0.08
	ū		0.02
		Initial condition (eq. (32))
	C a	abcissa of the center	0.25
	C γ σ 2	ordinate of the center variance	0.15 0.0025

Table 2 .

 2 Values of the local parameters (one per follicle) used in Figures 4 to 8

http://www.pas.rochester.edu/~bearclaw/

We use this index notation for the use of the forthcoming Multiresolution development

This simulation was performed for technical purposes and do not meet exactly the biological specifications. The precise numerical calibration of the model will be set up in future work.

Acknowledgment

The authors would like to thanks

• Pascal Joly, Philippe Parnaudeau and Kyrill Pichon Gostaff 4 for their help concerning parallel computing, • the CIRM 5 for the excellent working conditions, • the INRIA Large Scale Initiative Action REGATE 6 for the funding of the project.


Jean-Louis Combes
email: l.combes@u-clermont1.fr

Mireille Ntsama Etoundi
email: mireille.ntsama_etoundi@u-clermont1.fr

Thierry Yogo
email: yogout@gmail.com

Patrick Plane

Catherine Araujo Bonjean

Les auteurs

Keywords: Household consumption, food price shocks, vulnerability, aid, remittances

Les commentaires et analyses développés n'engagent que leurs auteurs qui restent seuls responsables des erreurs et insuffisances.

Introduction

The recent increase in food prices over previous years has raised concerns about the risk of the resurgence of undernourishment and hunger around the world, and mainly in poor countries. According to the World Bank (2011), global food prices remain high, partly due to increasing fuel prices as well as speculative behaviours in the commodity markets, and the World Bank's food price index has once again reached its peak from 2008. 1 Moreover, the World Bank study stressed that, since June 2010, an additional 44 million people have fallen below the $1.25 poverty line as a result of higher food prices. Simulations show that a further 10% increase in the food price index could lead to 10 million people falling into poverty, and a 30% increase could plunge 34 million people into poverty.

Notwithstanding the recognised damaging effects of food price shocks on poverty and welfare in developing food-importing countries, little is known about the role that international capital flows can play in dampening the effects of food price shocks. There are several reasons behind the interest in the effect of international capital flows, such as foreign aid and international migrant remittances. 2First, developing countries are among the top recipients of official development assistance (ODA) and remittances (World Bank Migration and Remittances Factbook, 2011) and there are existing studies on the countercyclical responsiveness of aid and remittances when countries are facing various types of shock, such as exchange rate shocks, natural disasters, conflict, oil shocks and financial crises [START_REF] Ratha | Workers' remittances: an important and stable source of external development finance, Remittances: Development Impact and Future Prospects[END_REF][START_REF] Yang | Coping with disaster: The impact of hurricanes on international financial flows, 1970-2002[END_REF][START_REF] Lueth | Are Workers' Remittances a Hedge Against Macroeconomic Shocks? The Case of Sri Lanka[END_REF][START_REF] Mohapatra | Remittances and natural disasters : ex-post response and contribution to ex-ante preparedness[END_REF][START_REF] Dabla-Norris | Business Cycle Fluctuations, Large Shocks, and Development Aid: New Evidence[END_REF][START_REF] David | How do international financial flows to developing countries respond to natural disasters?[END_REF]. 1 Several factors lie behind the recurrent rises in food prices. These drivers include severe weather events in key grain-exporting countries, the broad-based increase in agricultural commodity prices in 2010, which increased the competition for land and other inputs; and the link between higher oil prices and biofuels.

Second, there are very few papers which examine the role of international capital flows as shock absorbers in developing countries [START_REF] Collier | Aid, shocks, and growth[END_REF]Chami et al., 2009;[START_REF] Craigwell | Economic volatility and remittances[END_REF][START_REF] Combes | Remittances and household consumption instability in developing countries[END_REF]. These papers essentially focus on the contribution of international capital inflows to macroeconomic stabilisation and coping with various types of shock. [START_REF] Collier | Aid, shocks, and growth[END_REF] found that the adverse effects of negative export price shocks can be mitigated by broadly contemporaneous increases in aid. [START_REF] Combes | Remittances and household consumption instability in developing countries[END_REF] showed that remittance inflows help to reduce the destabilising effects of natural disasters, agricultural shocks, financial and banking crises, discretionary fiscal policies and exchange rate volatility on household consumption per capita. One principal limitation of these studies is that they have neglected to compare the effectiveness of remittances and foreign aid following a shock, and especially food price shocks in developing countries. This issue is important, since there is a fear that foreign aid will be crowded out by the large remittance flows observed nowadays. 3 Hence, it seems worthwhile to compare the effectiveness of these two flows in times of shock before making any judgments.

Third, there are very few macroeconomic papers which examine the effects of the rise in food prices by distinguishing between the potential effects on the most exposed countries (the vulnerable ones) and those in which there are not necessarily any negative effects on welfare (the less vulnerable ones). One exception is [START_REF] Kamgnia | Political Economy of Recent Global Food Price Shocks: Gainers, Losers and Compensatory Mechanism[END_REF], who addressed the effects of the level of food prices on various outcomes (undernourishment, agricultural production, current account and government spending), and differentiated the effects according to the level of vulnerability to food price shocks in Sub-Saharan African countries.

3 According to [START_REF] Grabel | Remittances: Political Economy and Developmental Implications[END_REF], skeptics regarding ODA and international aid bureaucracies have embraced remittances as part of what [START_REF] Adelman | The privatization of foreign aid: Reassessing national largess[END_REF] approvingly calls the new "privatized foreign aid" and what the Financial Times (cited in Adelman) terms the "diaspora that fuels development." In this view, remittances are superior to traditional (public) ODA, because they have little to no overhead, they are not subject to misuse by state officials and they efficiently and directly meet human needs in developing countries.

This paper measures the effects of food price shocks on both the level of household consumption per capita and the instability of the household consumption per capita growth rate in developing countries. In this vein, this paper explores the role of aid and remittance inflows in the mitigation of the effects of food price shocks in the recipient economies.

Using a large sample of developing countries observed over the period 1980-2009 and after mobilising dynamic panel data specifications, the econometric results yield three important findings. First, food price shocks significantly affect both the level and the instability of household consumption in highly vulnerable countries. Based on the previous work of de [START_REF] De Janvry | The global food crisis: Identification of the vulnerable and policy responses[END_REF], the vulnerability of countries to food price shocks was assessed by computing a continuous index which aggregates three main dimensions: the level of underdevelopment; high food dependency and a high food import burden. Second, remittance and aid inflows significantly dampen the effect of food price shocks in the most vulnerable countries. Third, a lower remittance-to-GDP ratio is required in order to fully absorb the effects of food price shocks compared to the corresponding aid-to-GDP ratio which is required.

The remainder of the paper is as follows. Section 2 computes the vulnerability index, the food price shocks and provides some stylised facts about the vulnerable countries. Section 3 presents the econometric models and discusses the main results. This paper concludes in Section 4 with policy implications.

Differential effects of food price shocks according to the level of vulnerability: Some stylised facts

Identifying vulnerable countries

According to de [START_REF] De Janvry | The global food crisis: Identification of the vulnerable and policy responses[END_REF], countries are vulnerable if they meet the following three criteria: (1) high food dependency; (2) a high food import burden and (3) low income.

The first criterion highlights the importance of food in the basket of goods consumed by the representative household in a given country. Hence, the larger the share of food items in the basket, the harder the household will be hit by an increase in food prices. We approximate food dependency by the share of total food imports in the total household consumption. The second criterion stresses the strong dependency of a country on the rest of the world, including for its supply of food items. Thus, net food-importing countries would be naturally more vulnerable to increases in food prices because such increases would cause their current account to deteriorate considerably. In this paper, we measure this variable by the ratio of food imports to total imports. Finally, the third criterion underlines the capacity of a country to constitute food safety nets for domestic consumers, should price shocks be prejudicial to them. We measure this capacity by the level of income as approximated by GDP per capita.

Using the criteria of de Janvry and Sadoulet (2008), we resorted to principal component analysis (PCA) in order to build our vulnerability index. Thus, the vulnerability index is a combination of the following variables: the ratio of food imports to total household consumption; the ratio of total food imports to total imports of goods and services and the inverse of the level of GDP per capita. The latter transformation is achieved in order to ensure that the level of development is negatively correlated to the degree of vulnerability to food price shocks. Finally, the vulnerability index is rescaled so that it ranges between 0 and 10, with higher values corresponding to high levels of vulnerability. The calculation of the vulnerability index is made for the period 1980-2009. The statistical summary of the PCA is presented in the appendices (see Tables A1 andA2 in Appendix A.).

Overall, the first two principal components explain more than 89% of the total variance and are significantly correlated with the three main variables used to build the vulnerability index.

Once the vulnerability index is built, it is useful to see how it is distributed among income groups and regions worldwide.

Figure 1 presents the distribution of vulnerability to price shocks among developing regions.

The figure shows that the vulnerability index is greater for Sub-Saharan African countries. This result is, however, intuitive. In fact, most of these countries are highly dependent on food imports and do not have enough income to build safety nets to cope with food price shocks. Note : The vulnerability to food price shocks index is the aggregation of three variables: the inverse of GDP per capita, the level of food imports as percentage of total imports of goods and the level of food imports as percentage of total household final consumption. The principal component analysis has been used as the technology of aggregation and the resulting vulnerability index has been rescaled to be between 0 and 10 with higher value indicating a strong level of vulnerability. Looking at income groups, it appears that low-income countries are the most vulnerable, as shown in Figure 2. This is due mainly to their inability to build a safety net, such as buffer stocks, as underlined above. Overall, the distribution of the vulnerability index among the developing regions and income groups shows that low-income countries are the most vulnerable. However, among this group, Sub-Saharan African countries are the least secure.

Computing the food price shock variable

In order to compute the food price shock variable, we follow the methodology developed by [START_REF] Deaton | International Commodity Prices, Macroeconomic Performance, and Politics in Sub-Saharan Africa[END_REF] and used by [START_REF] Dehn | Commodity price uncertainty and shocks: implications for economic growth[END_REF] and [START_REF] Collier | Aid, shocks, and growth[END_REF].

Specifically, we follow two steps: (1) the food price index is computed following the work of [START_REF] Deaton | International Commodity Prices, Macroeconomic Performance, and Politics in Sub-Saharan Africa[END_REF] and (2) food price shocks are identified using the basic forecast model of [START_REF] Collier | Aid, shocks, and growth[END_REF].

Let t i

FP , be the food price index in country i for the year t ; t j P , be the world price of food item j at time t and

• j i
w , be the country-specific weighting of the food item at the base year, so that we have the following geometrically-weighted structure of the commodity price index:

∏ = • = 6 1 , , , j t j t i j i w P FP where • j i
w , is the value of food item j in the total value of all commodities 6 n = for the constant base period. In this paper, the basket of goods is made up of six commodities: maize; milled rice; soybean oil; soybeans; refined sugar and wheat. These commodities are considered as the most commonly imported foods in the world (FAO, 2011). 4 The weighting item is obtained using the following formula:

∑ = • • • • • = 6 1 , , , j j i j j i j j i Q P Q P w where • , , j i Q
represents the import value of commodity j at the base year. The year 1995 was chosen as the year for which to construct country-specific commodity import weights, because it allows the inclusion of the former Soviet Union countries and maximises the number of observations of food imports disaggregated by the products in the sample [START_REF] Burke | Do output contractions trigger democratic change?[END_REF]. Global food prices by product are drawn from the International Monetary Fund World Economic Outlook database, and the data on food imports are drawn from FAO TradeStats online.

The second step consists of computing the shock variable. Following [START_REF] Collier | Aid, shocks, and growth[END_REF], shocks are located by differencing each country's aggregate real commodity price index series in order to make it stationary, removing predictable elements from the stationary process and normalising the residuals. The forecasting model used to identify shocks is estimated for each country separately, as follows:

t i t i t i t i FP FP t FP , 2 , 2 1 , 1 1 0 , ε θ θ α α + + ∆ + + = ∆ - -
with i denotes the country and t the year. The residuals from the equation above ( t i, ε

) ) are normalised by subtracting their mean value and dividing by their standard deviation. Food price shocks are the positive observations of the residuals t i, ε

) . In other words, the negative observations of the normalised residuals t i,

ε

) are replaced by zeroes in the database. More formally, the food price shock t i S , is written as follows:

] 0 [ 1 , , > = t i t i S ε )

Comparing the level of household consumption in times of shock between vulnerable and non-vulnerable countries

This subsection compares the effects of shocks on household consumption between vulnerable and non-vulnerable countries. For this purpose, the relative deviation of the real household consumption per capita during food price shocks compared to a situation without a food price shock is computed. High refers to a level of vulnerability which exceeds the 75 th percentile of the vulnerability variable, intermediate refers to a level of vulnerability between the median value and the 75 th percentile of the vulnerability variable, and low refers to a level of vulnerability below the median value of the vulnerability variable in the sample.

One may notice that the household consumption deviation value is negative for the most vulnerable countries, slightly positive for the intermediates and fairly high for the less vulnerable countries. In other words, countries which are highly vulnerable exhibit a decrease in consumption in times of shock. Those which are less vulnerable enjoy an increase in consumption relative to times when no shock has occurred. This is due mainly to the fact that most of these countries are net food exporters, and therefore benefit from the increase in prices. Figure 4 shows that in times of shock, the most vulnerable countries receive more ODA than the other groups of vulnerable countries. The increase in ODA relative to times of an absence of shocks stands at 13.6%, while this figure is 8.4% for the "intermediate" and 1.5% for the "low" group of countries. Levels of vulnerability to food price shocks Note : The foreign aid shock records the relative deviation (in percentage) of the aid-to-GDP ratio in the time of food price shocks compared to a situation without a food price shock in each country.

High: refers to a level of vulnerability exceeding the 75 th percentile of the vulnerability variable, Intermediate: refers to a level of vulnerability comprising between the median value and the 75 th percentile of the vulnerability variable, Lows: refers to a level of vulnerability below the median value of the vulnerability variable in the sample.

Source : FAO-Stats, World Bank Development Indicators and authors' construction.

In the same vein, Figure 5 presents the deviation in remittances relative to the absence of shocks. This figure suggests that in the event of shocks, the highly vulnerable countries once again receive more remittances than the "intermediate" and "low" group of countries.

However, what is striking in this case is the fact that the increase in the amount received is extremely high, standing at 143.9% relative to situation without shocks.

Do aid and remittances dampen the effect of food price shocks? Econometric models

The following models are specified in order to test the impact of food price shocks on household consumption.

a) Models of the effects of food price shocks on household consumption

The first models describe the effect of food price shocks on the level and instability of the real household consumption per capita.

τ τ τ τ τ τ ε η θ β ρ , , 1 , 1 , , i i i i i i u S c c + + + + ′ + = - X (1) τ τ τ τ τ τ ε η φ γ σ ϕ σ , , 1 , 1 , , i i i i i i u S + + + + ′ + = - Z (2)
with τ , i c and τ σ , i representing the level of real consumption per capita and the instability of the real household consumption per capita growth rate (both expressed in logarithmic terms), respectively. X and Z represent the vectors of the control variables. S represents the average number of positive food price shocks in each country during each period. There may be a concern that households anticipate food price shocks and increase their consumption (if households tend to build buffer stocks in order to cope with future shocks). This can lead to biased estimates of the effects of food price shocks on consumption. However, this is less likely to occur, due to the fact that we are using data averaged over several years. Indeed, this overshooting with regard to consumption would be easily observed with high-frequency data (especially monthly data). . In other words, food price shocks, on average, have no significant effect on the level of household consumption, except in a context of high vulnerability to food price shocks. In order to account for this heterogeneity, models 1 and 2 are modified to include a nonlinear effect. In formal terms, the following models are estimated:

τ τ τ τ τ τ τ τ τ ε η θ θ θ β ρ , 1 , 3 1 , , 2 , 1 , 1 , , * i i i i i i i i i u Vul Vul S S X c c + + + + + + ′ + = - - - (3) τ τ τ τ τ τ τ τ τ ε η φ φ φ γ φσ σ , 1 , 3 1 , , 2 , 1 , 1 , , * i i i i i i i i i u Vul Vul S S Z + + + + + + ′ + = - - - (4)
The hypotheses tested are:

0 1 ≥ θ , 0 2 < θ and 0 1 ≤ φ , 0 2 > φ
. Models 3 and 4 allow the computation of the threshold of vulnerability to food price shocks, beyond which the detrimental effects of the shocks on welfare can be observed. From model 3, the vulnerability threshold is given as:

2 1 * 1 , 2 1 , , 0 = θ θ θ θ τ τ τ - = ⇒ = + ∂ ∂ - Vul Vul S c i i i
From model 4, we can see:

2 1 * 1 , 2 1 , , 0 = φ φ φ φ σ τ τ τ - = ⇒ = + ∂ ∂ - Vul Vul S i i i
For the model of the level of household consumption (model 1), the basic set of control variables includes the lagged dependent variable, per capita income, the instability of the growth in per capita income, the age dependency ratio and the rural population. Positive signs are expected for the lagged dependent variable and per capita income. The instability of the per capita income is expected to decrease the level of consumption, owing to the accumulation of precautionary savings, the decrease in private investment and the unavailability of jobs. The age dependency ratio is also expected to lower the level of consumption per capita when active people take care of less active individuals. Finally, in rural areas which are characterised by a low level of financial development and higher poverty rates, the share of consumption in household budgets would be higher. One could therefore expect a positive correlation between the rural population and the level of household consumption.

The equation for the instability of household consumption per capita is similar to those proposed by [START_REF] Herrera | Public expenditure and consumption volatility[END_REF], [START_REF] Craigwell | Economic volatility and remittances[END_REF] and [START_REF] Combes | Remittances and household consumption instability in developing countries[END_REF]. The variables that are expected to be positively correlated with the instability of household consumption are: the instability of the GDP per capita growth rate; government size and trade openness. The private credit ratio, the level of economic development and financial openness are expected to be negatively correlated with the instability of consumption in developing countries.

b) Models of the stabilising effects of foreign aid and remittances

Several specifications are adopted in order to test the hypothesis that foreign aid and remittances act as food price shock absorbers in the recipient economies. If this hypothesis is proven to be true, one would observe a decreasing marginal effect of food price shocks on household consumption as remittance and aid inflows rise. This would be the case specifically in the most vulnerable countries. In order to test this hypothesis empirically, this paper proceeds in two steps. First, the stabilising role of aid and remittances is evaluated using the entire sample of countries. We expect a non-significant effect of aid and remittance inflows.

Next, the sample of countries above the previously computed threshold of vulnerability (Vul*) is used to identify the stabilising effects of aid and remittances during food price shocks. As previously noted, the outcome variables are the level and the instability of household consumption per capita. The same matrix of control variables is also retained, along with the identification strategy through system-GMM. In formal terms, this appears as:

τ τ τ τ τ τ τ τ τ ε η θ θ θ β ρ , , 6 , , 5 , 4 , 1 , , * i i i i i i i i i u R R S S X c c + + + + + + ′ + = - (5) τ τ τ τ τ τ τ τ τ ε η θ θ θ β ρ , , 9 , , 8 , 7 , 1 , , * i i i i i i i i i u A A S S X c c + + + + + + ′ + = - (6) τ τ τ τ τ τ τ τ τ τ τ τ ε η θ θ θ θ θ β ρ , , 14 , 13 , , 12 , , 11 , 10 , 1 , , * * i i i i i i i i i i i i u R A R S A S S X c c + + + + + + + + ′ + = - (7)
For the consumption instability model, the following equations are estimated:

τ τ τ τ τ τ τ τ τ ε η φ φ φ γ φσ σ , , 6 , , 5 , 4 , 1 , , * i i i i i i i i i u R R S S Z + + + + + + ′ + = - (8) τ τ τ τ τ τ τ τ τ ε η φ φ φ γ φσ σ , , 9 , , 8 , 7 , 1 , , * i i i i i i i i i u A A S S Z + + + + + + ′ + = - (9) τ τ τ τ τ τ τ τ τ τ τ τ ε η φ φ φ φ φ γ φσ σ , , 14 , 13 , , 12 , , 11 , 10 , 1 , , * * i i i i i i i i i i i i u R A R S A S S Z + + + + + + + + ′ + = - (10)
with R and A representing the remittance-to-GDP and aid-to-GDP ratios, respectively.

- = ⇒ = + ∂ ∂ R R S i i i 8 7 * , 8 7 , , 0 = θ θ θ θ τ τ τ - = ⇒ = + ∂ ∂ A A S c i i i 8 7 * , 8 7 , , 0 = φ φ φ φ σ τ τ τ - = ⇒ = + ∂ ∂ A A S i i i c

) Data and methodology

Data for the control variables are drawn from the World Development Indicators Online, except for the financial openness series, which are drawn from Chinn and Ito's publicly-available dataset. 6 The narrow definition of remittance is employed to record remittances. In other words, the remittance variable only records the money sent back by migrants who have been residing in the host country for at least one year. 7 Remittance data are normalised by country GDP series. The foreign aid variable records the total amount of the official development assistance (ODA) to developing countries as a percentage of GDP. 8 The estimation of the dynamic panel models 1 to 10 presented above with the use of the OLS estimator is inconsistent, since the lagged dependent variables are introduced alongside country fixed effects. This bias is a particular concern here owing to the short temporal dimension of the dataset used. The system-GMM estimator must therefore be implemented.

The equations in levels and the equations in first differences are combined in a system and estimated with an extended system-GMM estimator which allows for the use of lagged differences and lagged levels of the explanatory variables as instruments [START_REF] Blundell | Initial conditions and moment restrictions in dynamic panel data models[END_REF]. 9 The GMM estimations control for the endogeneity of some explanatory variables.10 6 The original series in [START_REF] Chinn | A new measure of financial openness[END_REF] dataset, named KAOPEN, contains both positive and negative values between -2 and 2. We added the value of 2 to all of the observations in order to get only positive values. 7 For a detailed discussion of the reasons why the narrower definition is suitable in empirical macroeconomic studies, see Chami et al. (2009). 8 The reader may wonder whether it would not be more appropriate to use disaggregated data on foreign aid, such as food aid or agricultural aid, instead of total aid. At least two reasons justify this choice. First, one objective of this paper is to compare the effects of remittances and foreign aid in the absorption of food price shocks in developing countries. Given that there are no disaggregated data on remittances according to their uses, resorting to total aid allows the comparison with remittances. Second, disaggregated data regarding foreign aid are relatively scarce and are only available for a few countries over a short period of time.

Two specification tests are used to check the validity of the instruments. The first is the standard Sargan/Hansen test of over-identifying restrictions. The second test examines the hypothesis that there is no second-order serial correlation in the first-differenced residuals. 11

Econometric results

This section begins by discussing the effect of food price shocks on the level and the instability of household consumption. Next, the analysis turns to the econometric results of the mitigating role played by aid and remittances.

a) Heterogeneity in the effect of food price shocks

Tables 1 and2 present the results of the effect of food price shocks on both the level and the instability of household consumption per capita. In each table, the first columns describe the impact of food price shocks on the corresponding outcome variable. Regardless of the dependent variable that is retained (the level or the instability of consumption), the hypothesis that food price shocks do not exert a statistically significant effect on consumption is not rejected (column 1 of Tables 1 and2).12 Column 2 of Tables 1 and2 presents the results of the nonlinear effect of food price shocks on the two outcomes, depending upon the level of vulnerability to food price shocks. The results suggest that food price shocks become statistically significant in explaining both the level and the instability of consumption once the models allow for interaction with the vulnerability index. In other words, the results reveal a marginal decreasing and significant effect of food price shocks on the level of household consumption per capita as the level of vulnerability 11 To deal with the well-known problem of instrument proliferation raised by the system-GMM estimator [START_REF] Roodman | A note on the theme of too many instruments[END_REF], the matrix of instruments is collapsed and the number of lags is always limited to a fix order.

increases. Regarding the model of consumption instability, the results highlight a marginal positive effect of food price shocks on the instability of consumption, which increases with the extent of vulnerability. From the results presented in column 2 of Tables 1 and2, one can compute the threshold of vulnerability to food price shocks beyond which their effect becomes critical for household consumption. The results of Tables 1 and2 indicate that the values of these thresholds do not differ between the two tables. Indeed, the threshold of the vulnerability index stands at 1.75 (Table 1) and1.64 (Table 2). The percentage of countries above this threshold stands between 53% and 58%. This suggests that countries that are located precisely above the median value of the vulnerability index are particularly concerned by the damaging consequences of food price shocks on consumption. For the rest of the countries in the sample (those which are nonvulnerable and especially those with a vulnerability index close to 0), food price shocks are likely to increase the level of household consumption, but not necessarily its instability (according to the value and the significance of the coefficient associated with the additive term of food price shocks in Tables 1 and2).

To sum up, the preliminary econometric investigations highlight the heterogeneity in the response of household consumption to food price shocks, which depends upon the extent of the country's vulnerability to food price shocks. The next task consists of examining the role played by foreign aid and remittance inflows in absorbing food price shocks.

b) The mitigating role of aid and remittance inflows

Tables 3 and4 present the results of the econometric specifications, allowing a nonlinear effect of food price shocks conditional on the level of foreign aid and remittance inflows. For each dependent variable, the econometric models are first estimated using the whole sample (columns 1, 3 and 5 of Table 3 andcolumns 1 and4 of Table 4). As expected, when the entire sample of countries (a mix of vulnerable and non-vulnerable countries) is used, aid and remittances have no statistically significant effect in dampening the impact of food price shocks, except in the model of the level of household consumption. However, this specification should be taken with a pinch of salt, given that the GMM specification tests are not satisfactory (column 1 of Table 3 and column 1 of Table 4).

The next task consists of estimating the models with the restricted sample of vulnerable countries, and especially the sample of countries with a vulnerability index above the median value of the variable.

The results in Table 3 suggest that when countries exhibit a degree of vulnerability above the sample median, remittance and foreign aid inflows have a strong dampening effect on the negative impact of food price shocks on the level of household consumption (columns 2, 4 and 6). Indeed, as expected, the coefficient of the additive term of food price shocks is negative and statistically significant, whereas the coefficients of the interaction between shocks and aid and remittance inflows are positive and significant. The results also highlight the fact that the remittance and aid-to-GDP ratios required for full absorption of the effects of food price shocks on the level of household consumption are 5% and 13%, respectively. For the remittance threshold, approximately 30% of countries are concerned. For the aid-to-GDP threshold, 24% of countries are concerned. This result reveals that a comparatively low ratio of remittances is needed in order to cope with food price shocks.

Regarding the model of consumption instability, the results of Table 4 suggest that when countries exhibit a vulnerability index above the sample median, foreign aid inflows have a powerful dampening effect on the negative impact of food price shocks on the instability of the household consumption growth rate (columns 5, 6, 7 and 8 of Table 4). For remittance inflows, countries must be highly vulnerable (above the 75 th percentile of the distribution of the vulnerability index) in order to identify a significant stabilising effect of remittance inflows (columns 3 and 8). In other words, the instability of household consumption is affected to a lesser extent by food price shocks thanks to remittance inflows in countries that The estimation method is the one-step System-GMM. Time effects are included in all the regressions. Robust T-statistics are below the coefficients. Instability is the 5-year standard deviation of the growth rate of the real household consumption per capita. Data are averaged over six nonoverlapping 5-year periods between 1980 and 2009. *** p<0.01, ** p<0.05, * p<0.1 are highly exposed to food price shocks. A tentative explanation for this result could be that countries differ in their use of remittances according to their location on the distribution of the vulnerability index. Remittances may be more likely to be devoted to the financing of consumption in countries that are highly dependent upon food imports compared to other countries. This argument is also supported by the fact that the destabilising effect of food price shocks increases significantly when the sample is restricted to the countries above the 75 th percentile of the vulnerability index. Indeed, the coefficient associated with the food price shock variable introduced additively increases from 0.37 to 0.82 between columns 2 and 3, and from 0.74 to 1 between columns 7 and 8 of Table 4.

The results also highlight that the remittance and aid-to-GDP ratios required for full absorption of the effects of food price shocks in countries located above the 75 th percentile of the vulnerability index are approximately 9% and 29%, respectively. For the remittance threshold, around 13% of countries are concerned. The percentage decreases to only 5% for the corresponding aid threshold. 13 As has been shown previously, a comparatively low ratio of remittances is in fact necessary in order to cope with food price shocks.

Concluding remarks

The dramatic increase in food prices over previous years has revived an increasing concern about food security, mainly for the most vulnerable countries. This paper examines the impact of food price shocks on household consumption according to each country's level of vulnerability to food price shocks. It also addresses the ability of foreign aid and remittances to mitigate the impact of food price shocks. 13 Indeed, Mozambique and Nicaragua are the sole countries for which the average aid-to-GDP ratio is above the threshold of 29% for the period 1980-2009. Based on a large sample of developing countries, observed over the period 1980-2009, two main results are derived. First, food price shocks significantly affect both the level and the instability of household consumption, especially in the most vulnerable countries. In the context of a high level of vulnerability, food price shocks reduce the level of real household consumption per capita, while fostering the instability of household consumption. Second, the results highlight that when countries exhibit a high degree of vulnerability, remittance and foreign aid inflows have a strong dampening effect on the impact of food price shocks on household consumption. Finally, the results suggest that a lower remittance-to-GDP ratio is required in order to fully absorb the damaging effects of food price shocks on consumption compared to what is required in terms of foreign aid. This paper has clear policy implications. The results outlined above suggest that remittances and aid should be increased in order to mitigate the dramatic effects of food price shocks in the most vulnerable countries. However, this could be a short-term policy. Indeed, in the long term, one should address the issue of vulnerability. One way in which to deal with this specific issue is to invest massively in agriculture, while increasing the level of diversification. Foreign aid and remittances could be at the heart of this process.

Figure 1 :

 1 Figure 1 : Vulnerability to food price shocks among developing regions (1980 -2009)

Figure 2 :

 2 Figure 2 : Vulnerability to food price shocks among income groups (1980 -2009)

Figure 3

 3 Figure 3 depicts the consumption deviation during food price shocks over the period 1980-

Figure 3 :

 3 Figure 3: Consumption deviations in time of food price shocks (1980 -2009)

Figure 4 :

 4 Figure 4: Official Development Assistance deviations in time of food price shocks (1980 -2009)

Figure 5 :

 5 Figure 5: Remittance deviations in time of food price shocks (1980 -2009)

iu

  represents the country fixed effects and τ η are the period dummies. i and τ are the country and the non-overlapping sub-periods spanning from 1980 to 2009, respectively,

  When models 5 to 10 are estimated with the restricted sample of countries exhibiting a vulnerability index above the critical threshold of Vul* (the most vulnerable countries in the sample), the following hypotheses hold: 6, 8 and 9, the thresholds of the aid and remittance-to-GDP ratios which allow the full absorption of food price shocks in the group of the most vulnerable countries are

  Note : The remittance shock records the relative deviation (in percentage) of the remittance-to-GDP ratio in the time of food price shocks compared to a situation without a food price shock in each country. High: refers to a level of vulnerability exceeding the 75 th percentile of the vulnerability variable, Intermediate: refers to a level of vulnerability comprising between the median value and the 75 th percentile of the vulnerability variable, Lows: refers to a level of vulnerability below the median value of the vulnerability variable in the sample.

		150
		143.9
	Remittance shocks	50 100
		0	8.8	4.7
		High	Intermediate	Low
		Levels of vulnerability to food price shocks

Source : FAO-Stats, World Bank Development Indicators and authors' construction.

Table 1 : Food Price Shocks Vulnerability and household consumption, GMM Dynamic Panel Data results

 1 The estimation method is the one-step System-GMM. Robust T-statistics are below the coefficients. Data are averaged over eight nonoverlapping 4-year periods between 1980 and 2009. *** p<0.01, ** p<0.05, * p<0.1.

	Dependent Variable: log real household consumption per capita	(1)	(2)
	Food price shocks	-0.000761	0.0201*
		(0.00452)	(0.0109)
	Food price shocks*Vulnerability index		-0.0114**
			(0.00567)
	Vulnerability index		0.0230
			(0.0320)
	lag dependent variable	0.704***	0.698***
		(0.155)	(0.205)
	log (GDP per capita)	0.370***	0.385**
		(0.110)	(0.179)
	Age dependency ratio	-0.00298**	-0.00311*
		(0.00133)	(0.00166)
	Rural population (%)	0.00333	0.00374*
		(0.00243)	(0.00197)
	GDP per capita growth volatility	-0.00830***	-0.00539
		(0.00271)	(0.00393)
	Foreign Aid-to-GDP	0.00937*	0.0108*
		(0.00550)	(0.00639)
	Remittances-to-GDP ratio	0.00635	0.00707*
		(0.00480)	(0.00363)
	Intercept	-0.579	-0.706
		(0.543)	(0.455)
	Observations	446	386
	Number of countries	82	81
	Joint test of Aid*food price shocks, P-value		0.13
	Vulnerability index threshold		1.75
	Percentage of countries concerned		53%
	Arellano-Bond Test of AR(1) P-value	0.03	0.05
	Arellano-Bond Test of AR(2) P-value	0.32	0.42
	Hansen OID test, P-value	0.2	0.1
	Number of Instruments	19	27
	Notes:

Table 2 : Food price shocks, vulnerability and the instability of household consumption

 2 The estimation method is the one-step System-GMM. Time effects are included in all the regressions. Robust T-statistics are below the coefficients. Instability is the 5-year standard deviation of the growth rate of the real household consumption per capita.

	Dependent variable:
	log of sd of real household consumption per capita growth	(1)	(2)
	Food price shocks	0.105	-0.156
		(1.429)	(-1.184)
	Food price shocks*Vulnerability		0.0954**
			(2.141)
	Vulnerability food price shocks		-0.129
			(-0.987)
	lag of dependent variable	0.208**	0.113
		(2.324)	(1.030)
	GDP per capita growth volatility	0.489***	0.453***
		(7.877)	(6.813)
	Government final consumption ratio	0.00258	0.0448*
		(0.119)	(1.727)
	Initial GDP per capita (log)	-0.0672	-0.0669
		(-1.327)	(-0.719)
	Initial private credit-to-GDP ratio	-0.00429** -0.00542**
		(-2.469)	(-2.236)
	Trade openness	0.00633***	0.00568*
		(3.800)	(1.759)
	Financial openness index	-0.322	-0.509
		(-1.126)	(-1.450)
	Financial openness ²	0.0489	0.0905
		(0.867)	(1.332)
	Intercept	0.892**	0.965
		(2.156)	(1.571)
	Observations	367	330
	Number of countries	90	89
	Joint significance of food price shocks coeff.,P-value		0.087
	Threshold level of the vulnerability index		1.64
	Number of countries above the threshold		52
	Percentage of countries above the threshold		58%
	AR(1), P-value	0.000	0.000
	AR(2), P-value	0.672	0.539
	Hansen OID, P-value	0.219	0.157
	Nb of instruments	22	26
	Notes: Data are averaged
	over six nonoverlapping 5-year periods between 1980 and 2009.
	*** p<0.01, ** p<0.05, * p<0.1

Table 3 : Remittances, Foreign Aid, Food prices shocks and household consumption, GMM Dynamic Panel Data results

 3

	Dependent Variable: log real household consumption per capita	Full sample	Vul>median Full sample Vul>median Full sample	Vul>median
		(1)	(2)	(3)	(4)	(5)	(6)
	Food price shocks	-0.0118	-0.0237**	-0.00401	-0.0395**	-0.00696	-0.0543**
		(0.00777)	(0.0110)	(0.00682)	(0.0183)	(0.0110)	(0.0239)
	Food price shocks*Remittance-to-GDP ratio	0.00495**	0.00531***			0.00229	0.00446*
		(0.00218)	(0.00178)			(0.00210)	(0.00251)
	Food price shocks*Aid-to-GDP ratio			0.000457	0.00304**	0.000300	0.00340**
				(0.000819) (0.00137)	(0.000898)	(0.00155)
	Remittance-to-GDP ratio	0.00258	-0.000716			0.00411	-0.000161
		(0.00349)	(0.00354)			(0.00338)	(0.00291)
	Aid-to-GDP ratio			-0.00117	-0.00507	0.00373	-0.00395
				(0.00341)	(0.00378)	(0.00503)	(0.00407)
	lag dependent variable	0.592***	0.491***	0.488***	0.521***	0.474***	0.495***
		(0.151)	(0.138)	(0.0995)	(0.129)	(0.115)	(0.143)
	log (GDP per capita)	0.440***	0.517***	0.555***	0.514***	0.528***	0.485***
		(0.100)	(0.126)	(0.0813)	(0.115)	(0.0916)	(0.141)
	Total population of 65 years old (%)	0.0132***	0.0176	0.00806**	0.0107	0.0128***	0.0192*
		(0.00327)	(0.0109)	(0.00378)	(0.0104)	(0.00362)	(0.0115)
	Rural population (%)	0.00300	0.00276	0.00353	0.00333*	0.00142	0.00146
		(0.00266)	(0.00284)	(0.00224)	(0.00200)	(0.00242)	(0.00265)
	GDP per capita growth volatility	-0.00624**	-0.00364	-0.00737**	-0.00435	-0.00553**	-0.00478
		(0.00255)	(0.00279)	(0.00316)	(0.00323)	(0.00240)	(0.00298)
	Intercept	-0.552	-0.414	-0.661	-0.536	-0.332	-0.116
		(0.586)	(0.652)	(0.460)	(0.414)	(0.540)	(0.608)

Table 3 (continued)

 3 The estimation method is the one-step System-GMM. Robust T-statistics are below the coefficients. Data are averaged over eight nonoverlapping 4-year periods between 1980 and 2009. *** p<0.01, ** p<0.05, * p<0.1.

	Full sample	Vul>median Full sample Vul>median Full sample	Vul>median

According to the World Bank's (2011) recent report, "Migration and Remittances Factbook," the total amount of remittances received by developing countries is three times the level of foreign aid, which positions remittances among the top external sources of finance in the developing world.

FAO (2011). Women in agriculture, closing the gender gap for development: The state of food and agriculture report.

For equation 1, data are averaged over eight four-year sub-periods. For equation 2, data are computed over six sub-periods of five years each.

This paper uses the one-step system-GMM estimator developed by[START_REF] Blundell | Initial conditions and moment restrictions in dynamic panel data models[END_REF] for dynamic panel data.

In all specifications, food price shocks, the lagged vulnerability index, period dummies, initial GDP per capita, initial private credit, financial openness, output growth volatility, the rural population ratio and the age dependency ratio are taken as strictly exogenous. The other variables are supposed endogenous.

The specification tests associated with the system-GMM specifications give satisfactory results and do not invalidate the dynamic panel specifications.

Iran, Islamic Rep. Zambia Kazakhstan Kenya Note: Vul refers to the country specific mean of vulnerability, Median is the median of the distribution of countries according to their level of vulnerability, and Vul 75 is the 75 percentile of the distribution of the vulnerability index.