Frédérique Clément

Frederique Clement

Multiscale mathematical modeling of the hypothalamo-pituitary-gonadal axis

Keywords: Multiscale, Mathematical models, ovarian follicle, GnRH surge, cell kinetics, hormonal rhythms

Although the fields of systems and integrative biology are in full expansion, few teams are involved worldwide into the study of reproductive function from the mathematical modeling viewpoint. This may be due to the fact that the reproductive function is not compulsory for individual organism survival, even if it is for species survival. Alternatively, the complexity of reproductive physiology may be discouraging. Indeed, the hypothalamo-pituitary-gonadal (HPG) axis involves not only several organs and tissues, but also intricate time (from the neuronal millisecond timescale to circannual rhythmicity) and space (from molecules to organs) scales. Yet, mathematical modeling, and especially multiscale modeling, can renew our approaches of the molecular, cellular and physiological processes underlying the control of reproductive functions. In turn, the remarkable dynamic features exhibited by the HPG axis raise intriguing and challenging questions to modelers and applied mathematicians. In this article, we draw a panoramic review of some mathematical models designed in the framework of the female HPG, with a special focus on the gonadal and central control of follicular development. On the gonadal side, the modeling of follicular development calls to the generic formalism of structured cell populations, that allows one to make mechanistic links between the control of cell fate (proliferation, differentiation or apoptosis) and that of the follicle fate (ovulation or degeneration) or to investigate how the functional interactions between the oocyte and its surrounding cells shape the follicle morphogenesis. On the central, mainly hypothalamic side, models based on dynamical systems with multiple timescales allow one to represent within a single framework both the pulsatile and surge patterns of the neurohormone GnRH (gonadotropinreleasing hormone). Beyond their interest in basic research investigations, mathematical models can also be at the source of useful tools to study the encoding and decoding of the (neuro-)hormonal signals at play within the HPG axis and detect complex, possibly hidden rhythms, in experimental time series.

Introduction

This article deals with the multiscale modeling and analysis of some dynamical processes arising within the hypothalamo-pituitary-gonadal (HPG) axis, with a special focus on the female reproductive axis. The HPG axis can be considered as the paragon of neuroendocrine axes, since it both concentrates all remarkable dynamics that can be exhibited by these axes and owns its unique specificities, as gonads are the only organs that host germ cells. It involves neuronal and non-neuronal cells spread across the hypothalamus (and connected with other brain areas including the cortex), the gonadotropic cells in the pituitary gland, and the gonads : ovaries in females, testes in males (see Figure 1). Via hypothalamic neurons, the reproductive function is subject to many environmental cues such as the daylength, food availability and social interactions, as well as to internal signals such as stress or metabolic status. If the conditions are favorable, specific hypothalamic neurons secrete in a pulsatile manner the gonadotropin-releasing hormone (GnRH), the "conductor" of the reproductive axis. The pulsatile GnRH secretion pattern ensues from the synchronization of the secretory activity of individual GnRH neurons. The release of GnRH into the pituitary portal blood induces the secretion of the luteinizing hormone (LH) and follicle-stimulating hormone (FSH) by the pituitary gland. The gonadotropins LH and FSH act on somatic cells within the gonads to support the development of germ cells and their endocrine gonadal function. In turn, hormones secreted by the gonads (steroid hormones such as androgens, progestagens and estrogens or peptide hormones such as inhibin) modulate the secretion of GnRH, LH and/or FSH within entangled feedback loops. In females, the GnRH secretion pattern dramatically alters once per ovarian cycle, in response to the time-varying levels of ovarian steroids, resulting in the GnRH surge characterized by massive continuous release of GnRH. The GnRH surge is responsible for ovulation, leading to the release of fertilizable oocytes from ovarian follicles. Also, the GnRH pulse frequency changes along the ovarian cycle, and exerts a differential control onto the secretion of gonadotropins. The modeling approaches that will be overviewed here are specially involved in the understanding of the triggering of the GnRH ovulatory surge by the hypothalamus and the development of ovarian follicles from initiation up to ovulation. Depending on the specific physiological issue addressed by the model, their starting point can be either middle-out (from the intermediary, mesoscopic level up and down to the other levels), bottom-up (from the microscopic to the macroscopic level) or top-down (from the macroscopic to the microscopic level). All these approaches have been introduced elsewhere and we refer the reader in particular to [START_REF] Clément | Multiscale modelling of follicular selection[END_REF][START_REF] Clément | Coupled somatic cell kinetics and germ cell growth : multiscale model-based insight on ovarian follicular development[END_REF][START_REF] Vidal | A dynamical model for the control of the GnRH neurosecretory system[END_REF][START_REF] Krupa | A network model of the periodic synchronization process in the dynamics of calcium concentration in GnRH neurons[END_REF] for a detailed exposition. Some additional materials on the model formulations are provided as on-line supplements ; they are yet not intended to replace the complete and rigorous presentations of the models. Similarly, we will not give any direct citation to the large corpus of bibliographic references on which the biological assumptions of the models were grounded ; we again refer to the original expositions of the models, where there are thoroughly commented [START_REF] Clément | Growth kinetics of the granulosa cell population in ovarian follicles : an approach by mathematical modelling[END_REF][START_REF] Echenim | Multi-scale modeling of the follicle selection process in the ovary[END_REF][START_REF] Clément | Foliation-based parameter tuning in a model of the GnRH pulse and surge generator[END_REF][START_REF] Clément | Mathematical modeling of the GnRH-pulse and surge generator[END_REF]. Also, even if the model formulations are based on generic biological principles that can be applied to many mammalian species, most of the numerical applications presented here have been undertaken in the ewe, which appears to be a very valuable animal model in our context, for several reasons that will be exposed throughout this article.

F. Clément

Multiscale modeling of the HHG axis 1

2 Multiscale modeling based on middle-out approaches : the instance of follicular development

Ovarian folliculogenesis is a unique instance of development still occurring during adulthood. It spans several months, starting from the time when primordial follicles (made up of an oocyte surrounded by a few flattened granulosa cells) leave the quiescent pool and initiate a process of growth and functional maturation ending up either by ovulation (release of a fertilizable oocyte), or (in most of the cases) degeneration at any stage of development. After initiation, follicular development can be separated into two distinct periods. During basal development, the morphological structure of the follicle settles progressively as an antral follicle (spheroidal structure with a central cavity and two tissular layers, the granulosa sheltering the oocyte within the cumulus oophorus and the theca delimitating the follicle from the ovarian cortex) that becomes more and more responsive to gonadotropins. During terminal development, the follicle is strictly dependent on gonadotropin supply and becomes a very efficient tissue for steroidogenesis. Large domestic species are particularly interesting to investigate the bases of follicular development.

Especially, in the ovine species, experimentalists can have access to a variety of data ranging from cell kinetics of granulosa cells or ultrasonographic in vivo monitoring of follicle growth, to endocrine time series of pituitary and ovarian hormones, and, in addition, there exist in several sheep strains natural mutations affecting the ovulation number and follicle physiology [START_REF] Scaramuzzi | Regulation of folliculogenesis and the determination of ovulation rate in ruminants[END_REF].

Deterministic and continuous spatio-temporal formalism for terminal development

The terminal part of follicular development corresponds to the latest stages where both the selection of ovulatory follicles and ovulation occur, and whose salient features can be summarized as :

-it is an hormonally-controlled process in which the basic functional events are the responses of granulosa cells to the pituitary hormones and especially FSH ; -the level of pituitary hormones is in turn tuned within a feedback loop involving the whole cohort of growing follicles, hence the whole population of granulosa cells amongst those follicles. Terminal development is thus an intrinsically multiscale process, where the granulosa cell is the pivot linking the lower, intracellular level, on which the signaling machinery operates to convert the hormonal signal into a cell's fate (progression along the cell cycle, differentiation or apoptosis), with the upper, tissular level, on which the coordinated evolution of cells is converted into a follicle's fate (ovulation or atresia). Moreover, at these stages, the morphodynamic changes are rather simple (they are almost limited to an increase in the antrum diameter), while the biochemical status of cells change steadily (through the expression of steroidogenic enzymes and LH receptors). Also, the number of follicular cells is great (on the order of hundreds of thousand or millions), so that the cell number can be considered as a real (rather than integer) value ruled by a continuous formalism. Putting these observations together, we chose to design a middle-out approach centered on the granulosa cell level and focused on the spatio-temporal evolution of the cell densities in follicles and their interactions (indirect coupling induced by the pituitary-ovarian feedback loop) [START_REF] Echenim | Multi-scale modeling of the follicle selection process in the ovary[END_REF][START_REF] Clément | Multiscale modelling of follicular selection[END_REF]. The cell density corresponds to the local repartition of cells on a functional space, where they are spread according to their position within or outside the cell cycle and level of terminal maturation.

F. Clément Multiscale modeling of the HHG axis

The resulting model is a specific instance of structured cell populations grounded on the following master equation, provided for any f th follicle of a N f -sized cohort :

∂φ i f ∂t + ∂g i f (u f )φ i f ∂a + ∂h i f (γ, u f )φ i f ∂γ = -λ i (γ, U )φ i f on Ω i , i ∈ {1, 2, 3} , f ∈ {1, . . . , N f }
The structuring (continuous) variables, a and γ, correspond respectively to the age and maturity of the follicular cells ; they define the 2D (functional) domain on which the cell population evolves. The i index corresponds the cell phases, which are delimited both horizontally and vertically by given ranges of age and maturity values (see Figure 2) ; Ω 1 , Ω 2 and Ω 3 are respectively associated with phase G1, phase SM (aggregating phases S,M and G2 of the cell cycle) and phase D (differentiated phase after cell cycle exit). The specific expression of the model terms may depend on the space variables either indirectly, through a phase-dependency (within a given phase, the value of the function is the same whatever the cell location), or directly, with an explicit dependency on the variable. In addition, these terms may depend on variables representing the endocrine control : U (t) (plasma FSH) and u f (t) (locally available FSH). φ i f (t; a, γ) represents the density of cells of age a and maturity γ in phase i at time t within follicle f . Within an elementary surface δaδγ, the density can be viewed as the local cell crowding. The transfer from one phase to another is governed by appropriate conditions defined on the internal boundaries of the domain and grounded on the continuity in cell fluxes, for instance a mitosis-induced doubling condition on the SM-G1 interface. Within the domain, cells are transported rightwards with a speed defined by the (time-varying) aging velocity g i f , and, either upwards or downwards according to the sign of the maturation velocity h i f . In contrast to the process of cell proliferation, that occurs punctually at the mitosis age and is thus embedded as a boundary condition, the cell death process is distributed over ages, so that it appears in the right hand side, where λ i is the apoptosis rate. More details on the model formulation are provided in the appendix (see Sketch of the model for terminal follicular development). The hormonal control is set dynamically and collectively from the feedback pressure exerted by the whole cohort of follicles, denoted by M (t). The contribution of each follicle M f (t) corresponds to the continuous equivalent of the weighted sum of its follicular cells over the whole domain (first-order moment of the density according to the maturity) :

M f (t) = γφ f dadγ, M (t) = N f f =1 M f (t)
The global control variable U (t) is a decreasing sigmoidal function of M (t), while the local control variable u f (t) is proportional to U (t), with a rate evolving as an increasing sigmoidal function of

M f (t).
The behavior of each follicle in response to this endocrine environment can be studied according to different macroscopic markers such as the total number of viable cells (see Figure 3), cumulative number of cells lost through apoptosis, global maturity (that can be interpreted as its steroidogenic ability), growth fraction (proportion of proliferative cells), and mitotic index. At the same time, the microscopic cell dynamics underlying these macroscopic outputs can be monitored in detail : repartition of cells into the different phases of the cell cycle, numbering of generations in term of cell divisions or age at cell cycle exit, heterogeneity amongst cells (see Figure 2), instantaneous cell loss. On the follicular level, the critical freedom degrees are the velocities of aging (progression along the cell cycle) and maturation, which determine both the need for and responsiveness to FSH and,

F. Clément

Multiscale modeling of the HHG axis 3 as a consequence, the proper use of the follicular proliferative resource (progressive switch between proliferation and differentiation ultimately setting the size of the pool of steroidogenically active cells) and the transit time through the vulnerability window where the follicle is the most sensitive to apoptosis-induced atresia.

Even if the focus of the model is on the ovarian side, it nevertheless accounts for the other levels of the HPG axis involved in follicle selection and ovulation, through the time-varying level of FSH (a direct function of the ovarian maturity M (t)) and triggering of the ovulatory surge (defined indirectly from a threshold on this maturity). Despite its simplicity, the multilevel character of the model, in addition to its multiscale character on the ovarian level, allows one to investigate the impact of the subtle endocrine interplay within the HPG axis, including the degree of pituitary sensitivity towards the inhibitory feedback of the ovarian estradiol and inhibin, and hypothalamic sensitivity towards the stimulatory feedback of estradiol.

Stochastic and discrete spatio-temporal formalism for basal development

After initiation, the basal development of ovarian follicles spans a long period where follicles are not strictly dependent on the gonadotropin supply, even if they may be gonadotropin-responsive before the transition to terminal development. Along these stages, the most striking changes concern the morphodynamics of follicles, with the progressive organization of the follicle structure that is settled first by the oocyte growth and increase in the number of granulosa cell layers, and then by the formation of distinct thecal layers and completion of the antrum, that leads to the distinction between the cumulus and mural granulosa cells.

Up to now, we have been interested in the "compact" phase of basal development (until the appearance of the first antral gaps), when paracrine interactions between the oocyte and its surrounding granulosa cells are prominent, and the shaping of the growing follicle is mainly due to the balance between the increase in the oocyte volume and the proliferation of granulosa cells leading to an increasing number of cell layers. Also, the starting number of granulosa cells is extremely low (on the order of ten), so that the cell fates have to be considered on both an individual (discrete) and stochastic ground.

Since the morphologic information are as relevant as the functional ones along this compact phase, we chose to design a model that would allow us to follow not only the balance between oocyte growth and granulosa proliferation, but also to trace the (radial and tangential) location of granulosa cells with respect to the time-varying boundary of the oocyte [START_REF] Clément | Coupled somatic cell kinetics and germ cell growth : multiscale model-based insight on ovarian follicular development[END_REF]. The model considers three interacting scales : (i) a microscopic, local scale corresponding to an individual cell embedded in its immediate environment, (ii) a mesoscopic, semi-local scale corresponding to either anatomical or functional subareas of follicles and (iii) a macroscopic, global scale corresponding to the whole follicle. These three scales are intricately merged on the dynamical ground, since the main events (cell division or displacement), as well as the oocyte growth law involve at least two different scales. The main freedom degrees are, on one side, the average cycle duration of granulosa cells within a given cell layer, that results from the influence exerted by diffusive proliferative factors emanating from the oocyte, and, on the other size, the sensitivity of the oocyte to trophic factors emanating from granulosa cells. The relative contribution of the oocyte and granulosa cell-derived factors is evolving constantly, since, not only the cell number increases, but also, as a consequence of oocyte growth, the volume of the layers, hence the maximal number of cells per layer. On the mathematical ground, the model is driven by a multiscale stochastic process, informing at the same time on (i) the (discrete) number of cells around the oocyte N t , that is incremented by F. Clément Multiscale modeling of the HHG axis one each time a cell division occurs ; (ii) the location (X k ) of each k th cell amongst the N t cells (k = 1, • • • , N t ), informing both on the radial distance from the oocyte surface (corresponding to the layer number) and the angular location within a cell layer, and (iii) the individual age (A k ) of each k th cell, assessed as the time elapsed since the division that gave birth to it (the age of both daughter cells is reset at division). The age variable evolves constantly since the cell gets older with time between two divisions, while the location changes passively due to the oocyte growth and actively whenever a displacement event occurs. The oocyte growth is represented by a differential equation with a deterministic part (intrinsic oocyte growth) and a stochastic part accounting for the trophic effect of the follicular cells, i.e. the weighted contributions of cells according to the cell layer they belong to (the smaller the layer number, the greater the contribution).

The law of evolution of the stochastic process underlying the model formulation is too tedious to be described here in detail ; complementary information are nevertheless provided in the appendix (see Sketch of the stochastic model for basal follicular development). It is sufficient to say that it operates as a counting process, that registers all the events of cell division and displacement affecting the population, so that both the cell number and cell repartition between layers change. More specifically, the cell number increases along the successive divisions, new layers are created and progressively filled with new cells or cells moving from the deeper layers. In the same time, the volume and capacity of the layers increase as the oocyte diameter increases. The timing and pattern of these changes are ruled by the expression chosen for the probability laws of cell division (an increasing function of cell age and decreasing function of distance) and cell displacement (where the probability of motion depends on the local tolerance to overcrowding and the level of crowding in the cell neighborhood). As a result, as illustrated in Figure 4, one can follow the model outputs on different scales : the temporal evolution of the cell number, oocyte diameter and follicular diameter (sum of the oocyte diameter and the cell layer depth) on the macroscopic scale, and the repartition within the follicle of either individual cells or a subpopulation of cells, respectively on the microscopic and mesoscopic scale.

A proper balance between oocyte growth and granulosa cell proliferation is required for normal morphogenesis in follicular development. In addition to reproducing the first stages of follicular development in wild-type situations [START_REF] Lundy | Populations of granulosa cells in small follicles of the sheep ovary[END_REF], the model also helps to explain situations of imbalance that may lead either to a greater than normal oocyte surrounded by fewer granulosa cells and layers, as observed in naturally occurring genetic mutations in sheep, or, on the contrary, a smaller oocyte trapped within a deep cell corona.

Multiscale modeling based on top-down/bottom-up approaches : the instance of the secretory pattern of GnRH

The ovarian cycle is driven by the finely tuned pattern of GnRH secretion, subject to the feedback exerted by estradiol and progesterone of ovarian origin. During each cycle, the pulsatile regime switches to the ovulatory surge leading to massive GnRH release [START_REF] Christian | The neurobiology of preovulatory and estradiol-induced gonadotropin-releasing hormone surges[END_REF], while, along the pulsatile regime, the pulse frequency is slowed down during the luteal phase and increases steadily until the surge triggering during the follicular phase. The generation of this complex secretory pattern starts on the level of individual GnRH neurons, whose electric and ionic activities are coupled with the downstream secretory activity. Individual neuron activities are coordinated on the level of neuron assemblies such as clusters, up to the whole population level. On each level, the neuronal activities F. Clément Multiscale modeling of the HHG axis are modulated by numerous afferences coming from regulatory neurons.

The cyclic transition from GnRH pulses to GnRH surge

The question of the effects of ovarian steroids (estradiol and progesterone) on GnRH neurons can be investigated from different angles according to the species. The most precise neuroanatomical studies have been performed in rodents (mice and rats). In contrast, physiological studies intending to dissect in time the effects of ovarian steroids have mostly been undertaken in domestic species, especially in the ewe. This species is particularly useful for studying GnRH secretion rhythms, since it has a large body size compatible with repeated sampling of pituitary portal blood and cerebrospinal fluid and further analysis of GnRH time series, the duration of its ovarian cycle (around 21 days) makes it easier to dissect the different steps in the temporal sequence of steroid action and, as for follicular development, it is closer to human ovarian physiology compared to rodents. Since we were primarily interested in the physiological impact of the GnRH secretory pattern (central control of ovulation), we tackled first the question of alternation between the pulsatile regime and the surge regime of secretion. We thus adopted a top-down approach, that was designed on the most macroscopic level, that of the populations of GnRH and regulatory neurons, keeping in mind that the representation of the populations could be refined later and separated into different subpopulations or functional clusters. Our approach is based on the interactions between an average GnRH neuron representing the whole population of secretory neurons, and an average regulatory neuron, representing the combined effect of the different populations of regulatory neurons [START_REF] Clément | Mathematical modeling of the GnRH-pulse and surge generator[END_REF]. The corresponding dynamical system reads :

εδ dx dt = -y + f (x) ε dy dt = a 0 x + a 1 y + a 2 + cX ε dX dt = -Y + g(X) dY dt = X + b 1 Y + b 2
where f (x) and g(x) are two cubic functions (hence they have two local extrema) :

f (x) = -x 3 +3λx and g(X) = -X 3 + 3µX.
Even if it is not as such multiscale in space, the model is nevertheless clearly multiscale in time, with 3 timescales O(1), O(ε), O(δε) (from the slowest to the fastest). Each subsystem (x, y) and (X, Y ) is indeed itself a slow-fast system with two timescales and in addition the (X, Y ) system is slower than the (x, y) one. Slow-fast systems are widely used in modeling for electrophysiology and neurosciences. They are well suited to representing dynamics characterized by sudden changes such as action potentials.

The faster system corresponds to the average activity of GnRH secreting neurons, while the slower one corresponds to the average activity of regulatory neurons. The x, X (fast) variables relate to neuron electrical activities, while the y, Y (slow) variables relate to ionic and secretory dynamics.

In each system, the fast and slow variables feedback on each other. The coupling between both systems is mediated through the unilateral influence of the slow regulatory interneurons onto the fast GnRH ones (cX term). The pulsatile release of GnRH is associated with the ionic dynamics through a thresholding effect (see appendix Sketch of the model for the GnRH secretory pattern F. Clément Multiscale modeling of the HHG axis along the ovarian cycle). With the appropriate choice of parameter values, the regulating system operates within an oscillatory dynamic regime accounting for the cyclic character of the steroid feedback exerted onto the hypothalamus along the ovarian cycle. The coupling term aggregates the global balance between inhibitory and stimulatory neuronal inputs onto the GnRH neurons. As a result, the secreting system alternates between an oscillatory regime, similar as that of the regulating system (yet with a much higher, as well as time-dependent frequency), where this system exhibits series of pulses, and a quasi-steady state regime corresponding to the surge mode. This model with multiple timescales offers a single dynamical framework for both the surge and pulse regime of GnRH secretion, and accounts for the qualitative (i.e. the right sequence of secretory events) and quantitative (i.e. the frequency, duration, amplitude of secretory events) specifications drawn from experimental studies, which amounts to embedding time-and dose-dependent steroid control within the model (see further explanations in [START_REF] Vidal | A dynamical model for the control of the GnRH neurosecretory system[END_REF][START_REF] Clément | Modeling the dynamics of gonadotropin-releasing hormone (GnRH) secretion in the course of an ovarian cycle[END_REF]).

The model is able to meet precise quantitative relations between the secretion signal features. Apart from the total duration of the ovarian cycle, which is expressed in physical time, these relations can all be expressed as ratios, regarding (i) surge duration over the whole cycle duration, (ii) the duration of the luteal phase over that of the follicular phase, (iii) pulse amplitude over surge amplitude, and (iv) pulse frequency in luteal phase compared to follicular phase (see Figure 5). The model can also be used to perform in silico experiments inspired from several experimental protocols, that all deal with the steroid control of GnRH secretion, yet in very distinct situations corresponding to different underlying neuroendocrine mechanisms : default of progesterone priming during the luteal phase, that affects the amplitude of the subsequent surge (luteal deficiency situation), surge blockade induced by administration of luteal levels of progesterone during the follicular phase, short-term effects of either progesterone or estradiol bolus administration on the pulse properties.

Synchronization in individual GnRH neurons

In the previous model, the synchronization of the ionic/secretory activities underlying the pulsatile character of GnRH secretion (apart from the surge) was taken for granted. The question of synchronization is as much a neuroscience issue than a reproductive one. It can be tackled from a bottom-up viewpoint to try deciphering the emergent network behavior from individual behaviors. In this new approach, what is taken for granted is the oscillatory character of the ionic (Calcium) dynamics of individual neurons, and the challenge consisted in reproducing events of global synchronization over a background of asynchronous activities [START_REF] Krupa | A network model of the periodic synchronization process in the dynamics of calcium concentration in GnRH neurons[END_REF]. Such events can be observed in ex vivo culture of olfactive placodes, that take advantage of the extra-cerebral embryonic origin of GnRH neurons and enable experimentalists to monitor electrophysiological markers such as action potential, intracellular calcium levels or GnRH release on a timespan of several hours. An intriguing feature of these placodes is their ability to exhibit synchronized, large amplitude calcium peaks that occur at the same frequency than the species-specific GnRH pulse frequency, as it has been well documented in rodents [START_REF] Terasawa | Intracellular Ca2+ oscillations in luteinizing hormone-releasing hormone neurons derived from the embryonic olfactory placode of the rhesus monkey[END_REF]. To account for both the mostly asynchronous individual dynamics and the recurrent, yet relatively rare, events of synchronization, we proposed the following network model :

             dx j /dt = τ -y j + 4x j -x 3 j -φ fall (Ca j ) dy j /dt = τ εk j (x j + a 1 y j + a 2 -η j φ syn (σ)) dCa j /dt = τ ε φ rise (x j ) - Ca j -Ca bas τ Ca        j = 1 . . . N dσ/dt = τ δεσ -γ(σ -σ 0 )φ σ 1 N N i=1 Ca i -Ca desyn

F. Clément Multiscale modeling of the HHG axis

where N is the number of neurons and φ rise (x j ), φ fall (Ca j ) are increasing sigmoidal functions (see also the appendix, Sketch of the model for calcium oscillations in embryonic GnRH neurons).

In the absence of the coupling term (η j φ syn (σ))), the first three equations of this system represent the dynamics of the electric activity x j , recovery variable y j and intracellular individual calcium Ca j in an individual neuron. The recipe for the system formulation is similar as that of the x, y subsystem in the previous section, with an additional variable Ca j interplaying with variable x j .

The common feature to every neuron is the excitable and rhythmic character of calcium dynamics, yet each neuron exhibits its own period (interpeak interval, IPI) and amplitude, within a range of admissible values. When the coupling term is activated, the mathematical structure of the model accounts for a phenomenon of synchronization between neurons, that, roughly speaking, is based on the forward neurons waiting for the others in the peak stage, so that all peaks occur within a short range of time and the average calcium value is higher than the individual peak amplitudes.

After each synchronized peak, there is a short silent episode and then the individual oscillations resume again in an asynchronous manner until the next event of synchronization (see Figure 6). The synchronization events occur on a much slower time scale (typically 50-60mn in rodents) than the individual peaks (typical IPI of 6-8 min) ; the time constant τ scales the time unit to physical time (min). The model is also able to reproduce additional and less frequent experimental observations, such as partial recruitment of cells within the synchronization process or the occurrence of doublets of synchronization, which appears to be sensitive to the distribution of the parameter tuning the sensitivity of a neuron to the network dynamics (η j in the second equation). Within this modeling framework, the coupling operating here is of volume transmission type ; the question remains open, on both the modeling and experimental grounds, of which neuronal connectivity could operate in a similar way as the global variable σ.

Further comments on the modeling of the HPG axis

The different organs involved in the HPG axis both process on their own sides and communicate with one another within entangled feedback loops. Hence, up to some extent, they can be considered separately in the most pertinent frameworks, provided that the scientific logic of neuroendocrine axes is kept in mind through the definition of endocrine inputs and outputs. For instance, the ovulation timing results from the coordination between two controlled processes acting on the ovarian level (follicle selection) and hypothalamic level (GnRH surge triggering). The former needs a middle-out modeling approach coped with specific partial differential equations (conservation laws), while the latter involves excitable neuronal dynamics, that can be handled by coupled nonlinear (relaxation) oscillators. Yet, in the long term, disposing of real-time interconnected, multiscale models of the different organic components of the HPG axis (gonadal function, hypothalamic and pituitary control) would be very useful tools to address physiological and clinical challenging questions and also to teach reproductive physiology interactively.

Up to now, we have tried to illustrate how multiscale models formulated from appropriate formalisms can help investigate some remarkable dynamic properties of the HPG axis in a relatively compact way. With such formalisms the links between different scales is explicit and mechanistic, which goes far beyond the statistical association between events or phenomena observed on different scales (as is often the case for molecular markers of diseases for instance). The other side of the coin is that the analysis and simulation of such models are not obvious and even raise open problems for the community of Applied Mathematics, that need to be solved before the models can be fully exploited

F. Clément Multiscale modeling of the HHG axis

on the biological ground (we refer the interested reader to [START_REF] Clément | Foliation-based parameter tuning in a model of the GnRH pulse and surge generator[END_REF][START_REF] Clément | Optimal control of cell mass and maturity in a model of follicular ovulation[END_REF][START_REF] Michel | Multiscale modeling of follicular ovulation as a mass and maturity dynamical system[END_REF][START_REF] Shang | Cauchy problem for multiscale conservation laws : Application to structured cell populations[END_REF]). This is especially true as far as the validation and quantitative calibration of the model (choice of the parameter values) are concerned, which are particularly challenging issues in a multiscale modeling approach, that still need methodological developments to be solved and even rigorously defined. Indeed, in each modeling approach presented so far, attention is paid not only to the ability of the model to reproduce already complex qualitative behaviors occurring simultaneously on different scales, but also to meet detailed specifications summarizing the available quantitative knowledge. Meeting those constraints cannot be only a matter of brute force optimization procedures ; it is also grounded on the results of the theoretical analysis of the model equations. Also, in contrast to the ambient context of big data management, the quantitative calibration of physiological multiscale models is confronted with the lack of quantitative and dynamic data [START_REF] Aymard | Cell-kinetics based calibration of a multiscale model of structured cell populations in ovarian follicles[END_REF] and the technical difficulties still hampering the temporal monitoring of reproductive systems in vivo. On the gonadal side, there is a need for retrieving accurate information about cell dynamics (duration of the cell cycle and its different phases, changes in cell kinetic indexes as mitotic index, labeling index, growth fraction . . .). On the central side, the equivalent of functional imaging in neuroscience is still missing in neuroendocrinology to allow experimentalists to follow neurosecretory events in a non invasive manner.

In the approaches described so far, we focused on the distal (gonadal) and central (hypothalamic) levels of the HPG. The intermediary (pituitary) level also raises interesting dynamical issues, that are worth being studied on their own, and also with the perspective of ultimately connecting models designed on different levels of the HPG. In males, where the dynamical issues are mainly limited to the control of LH pulse frequency by the interplay between GnRH and testosterone, a quite comprehensive understanding of the feedback loops relating the testes to the hypothalamo-pituitary complex is possible and amenable to modeling. In [START_REF] Keenan | A stochastic biomathematical model of the male reproductive hormone system[END_REF] for instance, a model based on stochastic differential equations proposes an integrative description of the synthesis and release of GnRH and LH under the control of testosterone. In females, the central and endocrine control of gonadotropin secretion is far more complex and results both in the subtile coordination between FSH and LH levels along the ovarian cycle and the GnRH-driven ovulatory surge. Drawing a mechanistic and physiological picture of processes in play would need to capture the basic dynamical principles underlying the decoding of GnRH signal features (frequency, amplitude, duration of pulses) by pituitary cells, such as (i) the differential control exerted by GnRH pulse frequency onto gonadotropin secretion, starting from the differential expression of the specific FSH or LH β-subunit (reviewed in [START_REF] Thompson | GnRH pulse frequency-dependent differential regulation of LH and FSH gene expression[END_REF]), (ii) the original signaling cascade of GnRH amongst G protein coupled receptors due to its lacking a C-terminal intra-cytoplasmic tail, (iii) the interplay between the endocrine control exerted by ovarian hormones and paracrine control exerted locally by activin, follistatin and inhibin, and (iv) the concurrent pronounced changes in GnRH receptor numbers. GnRH signaling in pituitary cells is largely studied both on the experimental and modeling grounds (see as significant instances among others [START_REF] Lim | Negative feedback governs gonadotrope frequency-decoding of gonadotropin releasing hormone pulse-frequency[END_REF][START_REF] Bertram | A mathematical model for the actions of activin, inhibin, and follistatin on pituitary gonadotrophs[END_REF]) and cannot be developed in the scope of this review. We just point here that, to our knowledge, none of the biochemically-designed models can distinguish the effect of an increase in the cumulative dose of GnRH from that of a genuine frequency increase ("pure frequency effect" obtained by compensating the frequency for duration and/or amplitude of GnRH pulses). This ascertainment motivated us to study, in a very simplified setup (feedforward signaling motifs), how a given dose of hormone can induce different outputs from the target system, depending on how this dose is distributed in time [START_REF] Fletcher | Interpreting frequency responses to dose-conserved pulsatile input signals in simple cell signaling motifs[END_REF] and we found that nonlinearity in the steady state input-output function of the system predicts the optimal input pattern. Understanding how such input-output functions can be an emergent property of realistic signaling networks remains a totally

F. Clément Multiscale modeling of the HHG axis

open question.

Beyond their purpose for basic research and mechanistic knowledge, dynamical models can also be used in the context of the HPG axis from the more practical viewpoint of data analysis, and especially for model-based analysis of time series. Experimental time series of hormonal levels in (neuro-)endocrinological studies are, most of the time, obtained from the peripheral blood system. Due to the sampling process and clearance from the blood, they may differ drastically from the instantaneous hormone release [START_REF] Clarke | Intensive direct cavernous sinus sampling identifies highfrequency, nearly random patterns of FSH secretion in ovariectomized ewes : combined appraisal by RIA and bioassay[END_REF]. Moreover, they are subject to the inherent noise introduced by the hormonal assay. Studies based on deconvolution tools (see e.g. [START_REF] Veldhuis | A review and appraisal of deconvolution methods to evaluate in vivo neuroendocrine secretory events[END_REF]) have attempted to reconstruct the theoretical signal from experimental time series without any possibility of biological validation. We have introduced another approach using dynamical models integrating the properties of secretion events and mimicking the experimental protocol to generate synthetic time series that reproduce the whole process leading from secretion to experimental time series, through sampling and hormonal assay. This approach allowed us both to explain and to validate our DynPeak algorithm for the dynamical analysis of luteinizing hormone (LH) rhythm [START_REF] Vidal | Dynpeak : An algorithm for pulse detection and frequency analysis in hormonal time series[END_REF], which is freely available at https://www-rocq.inria.fr/sisyphe/paloma/dynpeak.html. A similar approach also helped up to retrieve as much information as possible from long-term LH time series and to detect hidden (circannual) rhythms by comparing the time-frequency signatures of experimental and synthetic data [START_REF] Vidal | Endogenous circannual rhythm in LH secretion : insight from signal analysis coupled with mathematical modeling[END_REF]. The HPG axis involves three main organic levels : the hypothalamus, pituitary gland and ovaries, whose activities are coordinated by entangled endocrine feedback loops. During most of the ovarian cycle, the secretion pattern of the hypothalamic neuro-hormone GnRH (gonadotropin-releasing hormone) is pulsatile. The release of GnRH into the hypothalamo-pituitary portal blood induces the secretion from the pituitary gland of the luteinizing hormone (LH), that also follows a clear pulsatile pattern, and follicle-stimulating hormone (FSH). LH and FSH control the development of ovarian follicles and their secretory activity (as well as that of the corpus luteum). In turn, hormones released by the ovaries (steroid hormones such as progesterone and estradiol or peptide hormones such as inhibin) modulate the secretion of GnRH, LH and/or FSH. All along the ovarian cycle, GnRH pulse frequency adapts to the steroid environment (low frequency during the progesteronedominated luteal phase and high frequency during the estradiol-dominated follicular phase). Once per cycle, in response to the increasing levels of estradiol emanating from the pre-ovulatory follicles, the GnRH secretion pattern alters dramatically and a massive release occurs. The GnRH surge triggers in turn the LH surge that induces the ovulation of the selected follicles. (Adapted from [START_REF] Clément | Modeling the dynamics of gonadotropin-releasing hormone (GnRH) secretion in the course of an ovarian cycle[END_REF]). Top panel : Cell number in a single ovulatory follicle subject to biological specifications retrieved in the ovine species. The change in the cell number is represented as a function of time (bottom horizontal axis), and diameter (top horizontal axis), with initial time corresponding to a 1mm diameter. The left vertical axis is tipped with a unit of one million cells, while the right one marks the ratio of cell number increase with respect to the initial number. The red bullets correspond to the experimental observations, the black lines to interpolation curves, and the blue dashed lines correspond to simulated values. The experimental set combines different sources of data that were used to relate the cell number to the follicular diameter on one hand and the follicle diameter to time (or follicular age) on the other hand (details in [START_REF] Clément | Growth kinetics of the granulosa cell population in ovarian follicles : an approach by mathematical modelling[END_REF]). The insert shows the associated decrease in the growth fraction. Bottom panel : simulation of a cohort of follicles. In this instance, the cohort is made up of eight follicles, amongst which two ovulate with a different stabilized number of granulosa cells, while the others degenerate (Courtesy of Marie Postel).

F. Clément Multiscale modeling of the HHG axis

F. Clément

Multiscale modeling of the HHG axis 14 22) at initial time are the same in all cases, while the values of parameters κ1 and/or λ1 differ according to the genotype. ++ case: λ1 = 500, κ1 = 3.5 10 -5 (orange (in the electronic version) cross in Figure 10); BB case: λ1 = 500, κ1 = 1. 25 10 -4 (open pink (in the electronic version) square in Figure 10); II case: λ1 = 60.9, κ1 = 4. 4 10 -4 (open cyan (in the electronic version) diamond in Figure 10).

of κ 1 defining an interval [κ 1min , κ 1max ] such that we clearly got over-or undersized oocytes. Searching within this interval by dichotomy, we got another κ 1 value, such that the follicle belongs to neither the large oocyte domain nor the small oocyte one. More precisely, this means that, outside this interval, the normalized difference between the observed diameter, d sim O , and the expected (experimentally observed) di- The bottom panel is a sketch of the GnRH ouput from the pulse and surge generator model for an idealized species, that shows the qualitative and quantitative sequence of secretory events (low pulsatile regime in the luteal phase, high frequency regime in the follicular phase, triggering of the surge, resumption of pulsatility). The pattern shown here is close to what is observed in the rhesus monkey, with as long a follicular phase as a luteal one. The top model shows the GnRH output subject to specifications designed in the ewe species. The middle panel is a handmade sketch of the corresponding levels of progesterone and estradiol that feedbacks on GnRH secretion through a complex hypothalamic circuitry. (Adapted from [START_REF] Clément | Modeling the dynamics of gonadotropin-releasing hormone (GnRH) secretion in the course of an ovarian cycle[END_REF]).

F. Clément

Multiscale modeling of the HHG axis

F. Clément Multiscale modeling of the HHG axis

Sketch of the model for terminal follicular development

We provide here a more comprehensive formulation of the multiscale PDE model ; we nevertheless refer to the more mathematically-oriented articles cited in the main text for a rigorous presentation.

To alleviate the notations, we drop the i exponent initially used to make the phase-dependency explicit. Also, we rename the structuring variables a and γ by the classical space variables x and y.

In each follicle, the cell density is governed by the master equation

∂φ f ∂t + ∂(g f (x, y, u f (t))φ f ) ∂x + ∂(h f (x, y, u f (t))φ f ) ∂y = -λ(x, y, U (t))φ f ,
defined on the whole numerical domain (Ω 1 ∪Ω 2 ∪Ω 3 ) consisting of the successive G1 and SM phases of the cell cycles, and the single differentiated phase D. At initial time, the follicular cells are distributed according to φ f (0, x, y) = φ 0f (x, y) ; all cells are within the first cell cycle and they are desynchronized. Their maturity is uniformly distributed within a subpart of the range (0, γ s ).

The aging and maturation functions are formulated according to the following expressions (here and in the sequel, all parameters are real positive constants) :

g f (x, y, u f ) = γ 1 u f + γ 2 in phase G1 1 in phases SM and D h f (x, y, u) = τ h (-y 2 + (c 1 y + c 2 )(1 -exp( -u ū ))) in phases G1 and D 0 in phase SM
The rate of cell death is non negative only in a strip overlapping the top of phases G1 and bottom of phase D ([y - s , y + s ], hatched area in Figure 2), where it is defined by

λ(x, y, U ) = Λ exp -(y -y s ) 2 /ȳ × (U max -U ) /U max
This rate is maximal when both U (t) takes its lowest value U min (penalization due to poor hormonal supply), and y = y s (highest sensitivity to apoptosis at the cell cycle exit transition).

The N f equations in the PDE system {φ 1 , . . . , φ N f } are linked together through the control terms u f (t) and U (t), that depend on the first-order moments in maturity of the densities :

m 1 f (t) = yφ f (t, x, y)dxdy (follicular maturity) M (t) = N f 1 m 1 f (t) (ovarian maturity)
The plasma FSH level U (t) is defined by

U (t) = S(M (t)) with S(M ) = U min + U max -U min 1 + exp(c(M -M )) δ
The locally bioavailable FSH level u f (t) is defined proportionally to U (t) as

u f (t) = b f (m 1 f (t))U (t) with b f (m) = b 1 + 1 -b 1 1 + exp (-b 2 (m -b 3 )) F. Clément
Multiscale modeling of the HHG axis

The model is completed by appropriate boundary conditions on both the inner and outer boundaries.

The inner boundary conditions hold on the interfaces separating the whole domain into the different cell phases. They are defined as :

1. condition of cell flux continuity on the interfaces between phases G1 and SM (horizontal cell fluxes : φ f (t, x + , y) = g f φ f (t, x -, y)) and between phases G1 and D (vertical cell fluxes : φ f (t, x, y + ) = φ f (t, x, y -)) ;

2. condition of mitosis-induced doubling of the cell fluxes on the interfaces between phases SM and G1 (birth of two daughter cells from one mother cell : g f φ f (t, x + , y) = 2φ f (t, x -, y)) ;

3. waterproof condition between phases SM and D (φ f (t, x, y + ) = 0).

The formulation of the boundary conditions on the outer boundaries are based on the facts that (i) there is no cell influx from outside, (ii) the maximal maturity is bounded and (iii) the number of successive cell cycles is adapted to the timespan of the biological process (alternatively, boundary conditions as those defined on the SM/G1 interface can be applied as periodic conditions on the vertical outer boundaries).

F. Clément Multiscale modeling of the HHG axis

Sketch of the stochastic model for basal follicular development

D O (t) = D O (0) + t 0 F det (D O (s)) i w i N i(t)ds
The growth of the oocyte diameter D 0 (t) is ruled by a combination of a deterministic (red) part F det (intrinsic growth with saturation at a maximal diameter) and a stochastic (blue) part driven by the weighted contribution of granulosa cells (modulation of the growth slope) : all cells within a same layer, N i (t), have the same weight w i , which decreases as i, the layer number, increases. Note that the maximal number of layers is time dependent and that the exact expression of the equation has been simplified for the sake of clarity.

2. Direct effect of the oocyte on granulosa cells : probability of cell division

p div (t) = 1 -exp -A k (t)/λ i
The instantaneous division probability of the k th cell (amongst all N t (t) cells) is ruled by a combination of cell age A k (the older the cell, the greater the probability) and the proliferative effect of the oocyte translated into a cell-layer dependent value of the average cell cycle duration (parameter λ i increases with the layer number i).

F. Clément Multiscale modeling of the HHG axis 20

Interaction-based shaping and growth of the follicle

As a result of the cell division events, the number of cells N t (t) increases in a incremental way (one cell more each time a cell division occurs). On the other side, the oocyte diameter increases and the capacity of each layer increases. The balance between both processes sets the level of crowding in the cell neighborhood and rules the probability for cell displacement :

p disp (L (i,j) t , t) = 1 1 + e - d i,j -µ σ L (i,j) t
represents an elementary volume element within the follicle (such as the red-filled areas in the center panel below, with i the cell layer and j the angular location within the layer). The t index recalls that this volume inflates with oocyte growth. d i,j is the local cell density within this element (number of cells per unit volume), µ is a reference crowding level (in the case when cells are assimilated to incompressible spheres fulling the whole available volume), while σ is a parameter of tolerance to overcrowding. Cell motion is locally isotropic, yet with a short spatial range, so that, as the cell number increases, new layers are progressively filled and the layer number increases. Putting all these terms together allows us to describe the law of evolution of the cell population, by a stochastic point process, denoted by Z t , describing both the current total cell number N t , and, for each k th cell, the time elapsed since the last cell division (A k ) and its location with respect to the oocyte (X k ) :

Z t = Nt k=1 δ (X k (t),A k (t)) ,
where δ stands the Dirac delta function (a distribution function which is concentrated at a single point in space). At initial time, Z 0 = N 0 k=1 δ (X k (0),A k (0)) corresponds to the primary follicle stage with one single layer of granulosa cells surrounding the oocyte.

F. Clément

Multiscale modeling of the HHG axis

Sketch of the model for the GnRH secretory pattern along the ovarian cycle

The model consists of two coupled systems. Each system is described by the well-known FitzHugh-Nagumo model, a simplified version of the biophysical Hodgkin-Huxley model, initially designed to explain the ionic mechanisms underlying action potentials in the squid giant axon. The systems are formulated in a similar way except that one is faster than the other. The coupling is unilateral ; the system corresponding to an average GnRH neuron is forced by that corresponding to an average regulatory neuron.

                                       ε dX dt = -Y + g(X) dY dt = X + b 1 Y + b 2          Regulating system εδ dx dt = -y + f (x) ε dy dt = a 0 x + a 1 y + a 2 + cX y out = y(t)χ {y(t)>y th }                 
Secreting system This yields the global model :

" dx dt = y + f (x), " dy dt = a 0 x + a 1 y + a 2 + cX " dX dt = Y + g(X), dY dt = X + b 1 Y + b 2 ,
y out (t) = y(t) {y(t)>y th } .

An appropriate choice of c and a 2 values ensures that the output y out of this m (Figure 1.5) driven by the periodic oscillation of X along the relaxation lim Typical oscillatory pattern possibly exhibited by a system g(X) and f (x) are cubic functions ; typical expression can be f (x) = -x 3 + 4x.

All parameters are positive constants. Time separation constants ε and δ are "small" with respect to the other parameters (with a typical value close to 0.01) and define the different timescales of the model. These parameters also control some quantitative specifications. ε affects the pulse frequency ; it is on the order of the inverse of the number of pulses along a whole ovarian cycle and affects the relative duration of the pulsatile versus surge regime. δ controls the relative duration of a pulse with respect to the average interpulse interval. The parameters b 1 and b 2 are chosen so that (X, Y ) operates in an oscillatory regime ; (X, Y ) follows a so-called relaxation cycle alternating slow and fast parts. The slow parts correspond to the point (X, Y ) moving close to either the left or right branch of the graph of the function Y = g(X) in the (X, Y ) phase plane (teal and red parts marked by one arrow-head in the Figure above), while the fast parts correspond to almost horizontal jumps back and forth from one branch to the other (orange parts marked by two arrow-heads). In turn, the behavior of the secreting system is driven by the periodic behavior of the regulating system. The secreting system also displays an oscillatory (pulsatile) pattern when the coupling term cX is such that the intersection point between the cubic function y = f (x) and the straight line y = a 0 x + a 1 y + a 2 + cX remains on the middle branch (see Figure below). The shape of the relaxation cycle is similar to that of the regulating system, yet its period is much shorter.

F. Clément

Multiscale modeling of the HHG axis 22

The pulsatile release of GnRH is associated with the ionic dynamics through a thresholding effect (χ is the indicator function, its value is 1 when its argument is true and 0 if it is false) ; y th is the critical intracellular calcium concentration needed to trigger a pulse (excitation-secretion coupling). When this point passes through the local extrema of the cubic function, the system undergoes a dramatic dynamic change (a Hopf bifurcation) and it switches to the surge mode, that corresponds dynamically to a quasi-steady state (the (x, y) point tracks a steady-state whose coordinates vary with time). 1. As long as σ < σ on , each cell follows it owns rhythm ; the cells are asynchronous and the mean calcium level remains low. When σ exceeds σ 0 , φ syn (σ) is activated and the cells for which η j is large enough produce, all almost at the same time, a higher calcium peak than in the asynchronous period of the oscillation. The mean calcium level then exceeds Ca desyn which resets σ to a value close to σ 0 . (Adapted from [START_REF] Krupa | A network model of the periodic synchronization process in the dynamics of calcium concentration in GnRH neurons[END_REF]).

F. Clément

Multiscale modeling of the HHG axis 25

Figure 1 :

 1 Figure 1 : Overview of the female hypothalamo-pituitary-gonadal (HPG) axis and endocrine control of the ovarian cycle

Figure 2 :Figure 3 :

 23 Figure 2 : Sketch of the simulation domain and microscopic outputs of the multiscale model for terminal follicular development Top panel : Functional domain in age (abscissa) and maturity (ordinate). Two successive cycles are represented in the lower part of the domain ; the upper part corresponds to the differentiated area after cell cycle exit from phase G1 ; the hatched area delimitates the zone of the domain where cell loss can occur through apoptosis. Bottom panel : instance of cell repartition along a simulation, with the color code indicating the local amplitude of the density. On the selected snapshot, the cell density is distributed over two consecutive cell cycles and the passage through the SM-G1 interface has resulted in a mitosis-induced doubling of the density. The part of the density distributed in phase D corresponds to cells that have exited the cell cycles during the G1 phases of both the current and former cell cycles. (Courtesy of Marie Postel).

Fig. 6 .

 6 Fig. 6. Simulation results corresponding to different oocyte and gran sheep follicles. Panels A, D, and G illustrate the relationships between fol panels B, E, and H illustrate the relationships between granulosa cell n and oocyte diameter; panels C, F, and I illustrate the relationships bet layers surrounding the oocyte and the number of granulosa cells. Pane of one simulation fitting with the data set of wild-type (++) sheep. Th B correspond to the pooled data sets ++ 1, ++ 2, and ++ 3, represent Figure 5. Panels D, E, and F: Results of two simulations corresponding low granulosa cell numbers (dotted curves) and fitting with the data set o oocytes and very low granulosa cell numbers (solid curves) and fitting with The points in panels D and E correspond to the data sets BB and II repre of Figure 5. Panels G, H, and I: Results of one simulation correspondin oocytes and high granulosa cell numbers. The oocyte diameter (33 µm) initial time are the same in all cases, while the values of parameters κ1 a to the genotype. ++ case: λ1 = 500, κ1 = 3.5 10 -5 (orange (in the ele Figure 10); BB case: λ1 = 500, κ1 = 1.25 10 -4 (open pink (in the elec Figure 10); II case: λ1 = 60.9, κ1 = 4.4 10 -4 (open cyan (in the electro Figure 10).

of κ 1 Fig. 6 .

 16 Fig. 6. Simulation results corresponding to different oocyte and granulosa cell growth rates in sheep follicles. Panels A, D, and G illustrate the relationships between follicle and oocyte diameter; panels B, E, and H illustrate the relationships between granulosa cell number, follicle diameter, and oocyte diameter; panels C, F, and I illustrate the relationships between the number of cell layers surrounding the oocyte and the number of granulosa cells. Panels A, B, and C: Results of one simulation fitting with the data set of wild-type (++) sheep. The points in panels A and B correspond to the pooled data sets ++ 1, ++ 2, and ++ 3, represented in panels A and B of Figure 5. Panels D, E, and F: Results of two simulations corresponding to (1) large oocytes and low granulosa cell numbers (dotted curves) and fitting with the data set of BB sheep (2) very large oocytes and very low granulosa cell numbers (solid curves) and fitting with the data set of II sheep. The points in panels D and E correspond to the data sets BB and II represented in panels C and D of Figure 5. Panels G, H, and I: Results of one simulation corresponding to an instance of small oocytes and high granulosa cell numbers. The oocyte diameter (33 µm) and cell number (22) at initial time are the same in all cases, while the values of parameters κ1 and/or λ1 differ according to the genotype. ++ case: λ1 = 500, κ1 = 3.5 10 -5 (orange (in the electronic version) cross in Figure 10); BB case: λ1 = 500, κ1 = 1.25 10 -4 (open pink (in the electronic version) square in Figure 10); II case: λ1 = 60.9, κ1 = 4.4 10 -4 (open cyan (in the electronic version) diamond in Figure 10).

Figure 4 :Figure 5 :

 45 Figure 4 : Multiscale outputs of the model for basal follicular development Leftmost panel : 3D-like views of microscopic outputs ; to observe the oocyte (big yellow cell), the cell mass constituted by granulosa cells (small green cells) of the follicle has been partially recessed. Left center panel : histological slices of follicles at different stages of development corresponding to model outputs in both flanking panels (courtesy of Danielle Monniaux). Right center panel : mesoscopic outputs showing the distribution of cells emanating from the same ancestor cell ; the intensity of blue staining is proportional to the proportion of clonal cells with respect to the whole population. Rightmost panel : Fitting of macroscopic outputs (A : oocyte diameter versus follicular diameter, B : follicle diameter versus cell number) to experimental data obtained in wild-type ewes [10]. Adapted from [2].

Figure 6 :

 6 Figure 6 : Calcium output in the GnRH network reproducing the experimental results observed in olfactive placodes [13]. Top panel : asynchronous regime in the absence of coupling. The dynamics of intracellular calcium is shown in three different neurons with slightly different periods and amplitudes (Courtesy of Alexandre Vidal). Bottom panel : recurrent synchronization of high amplitude peaks in the presence of coupling. Three synchronized peaks occur within the simulation, amongst a network of fifty neurons (the same color can be used for several individual neurons). (Adapted from [4]).

1 .

 1 a P e ll u c id a GRANULOSA KITL : Kit-Ligand, BMP : Bone Morphogenetic Protein, GDF : Growth differentiation factor, GJA : gap-junction protein α Effect of the granulosa cells on the oocyte growth

5 g

 5 the oocyte => PDE model > SPACE/AGE model kground 2D view of the possible displacement for cell lying in a given (yellow) elementary volume of the follicle ; a schematic view of the layers is superimposed on an histological slice

Figure 1 . 4

 14 Figure 1.4 Limit cycle of the Regulating System and corresponding time trace of v and b2 are chosen so that the Regulating System admits an unstable singular point l X-nullcline surrounded by a limit cycle (left panel). Due to the slow-fast property of t of Relaxation type and the periodic X time trace is characterized by fast transitions be from Xmin to and slow decrease from Xmax to .

Figure 1 . 3

 13 Figure 1.3 Phase portrait of the GnRH Secreting System and output yout signal in the pulse regime (panels A-B) and surge regime (panels C-D-E). In each panel A and B, the cubic curve and the straight line represent the x-and y-nullcline respectively. The orange point corresponding to their intersection represents a singular point of the GnRH Secreting System. For I in an interval of values defined by the other parameters, the GnRH Secreting System admits an unstable singular point surrounded by an attractive limit cycle of relaxation type (panel A). Along this orbit, the generated y signal (grey pattern in panel B) is serrated and the thresholded yout signal (blue pattern in panel B) is pulsatile. For greater values of I, a stable singular point lies on the left branch of the cubic x-nullcline. Due to the slow-fast property of the GnRH Secreting System, the orbits first reach the x-nullcline quickly, and then follow it while tending to the singular point (panel C). Panels D and E show the time trace of y along the lower and upper orbits displayed in panel C respectively.

Figure 1 . 3 6

 136 Figure 1.3 Phase portrait of the GnRH Secreting System and output yout signal in the pulse regime regime (panels C-D-E). In each panel A and B, the cubic curve and the straight line represent the x-andThe orange point corresponding to their intersection represents a singular point of the GnRH Secret interval of values defined by the other parameters, the GnRH Secreting System admits an unstable si by an attractive limit cycle of relaxation type (panel A). Along this orbit, the generated y signal (gre serrated and the thresholded yout signal (blue pattern in panel B) is pulsatile. For greater values of I lies on the left branch of the cubic x-nullcline. Due to the slow-fast property of the GnRH Secreting reach the x-nullcline quickly, and then follow it while tending to the singular point (panel C). Panels trace of y along the lower and upper orbits displayed in panel C respectively.

Sketch of the model for calcium oscillations in embryonic GnRH neurons

Individual calcium dynamics

Calcium-induced hyperpolarization of the neuron membrane

Calcium dynamics induced by the polarization of the neuron membrane

τ Ca modulates the decay rate of calcium level to the baseline Ca bas

Embedding of the individual dynamics within a network of N neurons

σ is a global variable representing the network state acting on each neuron through a thresholded effect

and subject to a reset mechanism induced by each episode of calcium synchronization φ σ (u) = 1 1 + exp(-ρ σ u)

F. Clément

Multiscale modeling of the HHG axis


M El

H El Maroufy

Moulay Sultan

Slimane

Inference for Nonlinear Stochastic Equations with Random Effects. Application to Epidemic Model

We consider the maximum likelihood estimator of the parameters for stochastic differential equations with a drift that depends on random effects. An explicit formulae was given and the asymptotic properties are investigated when the random effects have Gaussian distribution. In this article, we deal with a nonlinear model, by using the approximate Log-likelihood instead of the exact one. We prove the asymptotic properties using existing relevant theorems when the drift term does not depend linearly on the random effects. An example of application in epidemiology is treated and discussed in detail.

Introduction

In biomedical research, studies in which repeated measurements are taken on series of individual or experimental animals play an important role. Models icluding random effects to model this kind of data enjoy an increasing popularity. These models account for both intra-individual variability, that is, occuring within the dynamics of each individual over time, and the variability occuring between subjects. The mixed effects SDEs have been introduced as an alternative to classical ODE-based models (Ditlevsen and De Gaetano, 2005b [4]; Overgaard et al., 2005 [11]; Donnet and Samson, 2008 [5]) in Pharmacokinetic/Pharmacodynamic modelling. SDEs with random effects have also proposed for neuronal data (Picchini et al., 2008 [13]).

Maximum likelihood estimation (MLE) of parameters of random effects, also called population parameters, is generally not straightforward as the likelihood or log likelihood function can rarely be expressed in closed form because of complexity of logarithm and integration. Approximations of the likelihood have been proposed, based on linearization [START_REF] Beal | Estimating population kinetics[END_REF] or Laplce's approximation (Wolfinger, 1993 [16]). Alternative methods have also been developed such as the SAEM algorithm (Kuhn and Lavielle, 2004 [9]). Maximum likelihood estimation in SDEs with random has been tackled in few papers. Ditlevsen and De Gaetano (2005a) show that in specific case of mixed effects Brownian motion with drift, the likelihood function can be explicitly derived, leading to explicit parameters estimators. For general mixed SDEs approximations of the likelihood have been proposed (Picchini et al., 2010[12]; Picchini and Ditlevsen, 2011 [START_REF] Picchini | Practicle estimation of a high dimensional stochastic differential mixed-effects models[END_REF]).

2 Model and assumptions [START_REF] Delattre | Maximum likelihood estimation for stochastic differential equations with random effects[END_REF] consider classical inference in the context of mixed-effects stochastic differential equations (SDE's) having the form

dX i (t) = b(X i (t), φ i )dt + σ (X i (t), φ i )dW i (t), (1)
where, for i = 1, ..., N , X i (0) = x i , is the initial value of the stochastic process X i (t), which is assumed to be continuously observed on the time interval [0, T i ],

where T i > 0 assumed to be known. Each X i models an individual with dynamics given by (1). W 1 , ..., W N are N independent Wiener processes, φ 1 , ..., φ N are R d -valued random variables which represent random effects. (φ 1 , ..., φ N ) are independent of Brownian motions and identically distributed with common distribution g(ϕ, θ)dν(ϕ). Here g(ϕ, θ) is density with respect (w.r.t for short) to a dominating measure ν on R r , for all θ, where θ ∈ Θ ⊂ R d is the unknown parameter of interest, which is to be estimated. Delattre et al., 2012 [2] impose regularity conditions that ensure existence of solutions of (1). Our interest is the estimation of the population parameter θ, where b(., .) and σ (., .) are real valued functions, but b(x, φ i ) is not necessary linear in φ i as considered in Delattre et al., (2012) [START_REF] Delattre | Maximum likelihood estimation for stochastic differential equations with random effects[END_REF]. Thus we construct the approximate likelihood function. To simplify notations, we write dϕ instead of dν(ϕ). We also consider firstly the one dimensional case where ϕ, x ∈ R, which is be extended in multidimensional case. Notationally,

P -→, D

=⇒ and

Pa.s -→, denote, respectively, convergence in probability, in distribution and almost surely under P. The following assumptions and notations will be needed throughout the paper.

Assumption (1) ( Smoothness of the coefficients )

(i) b(x, ϕ) ∈ C 1 (R 2 ,
R) and its x-derivative is bounded uniformly in ϕ.

(ii) σ (x, ϕ) ∈ C 2 (R 2 , R) and all its x-derivatives up to order 2 are bounded uniformly in ϕ.

Assumption (2) ( Non degeneracy of the diffusion )

There exists σ 0 > 0 such that, σ 2 (x, ϕ) ≥ σ 2 0 > 0, for all x, ϕ ∈ R

Assumption (3)

There exits L > 0 such that, b(x, ϕ) σ (x, ϕ)

+ |σ (x, ϕ)∂ xx σ (x, ϕ)| ≤ L for all x, ϕ ∈ R,
where ∂ xx f (x, ϕ) denotes the second derivative of f with respect to x.

Assumption (4) ( Law of random effects ) g(ϕ, θ)dϕ = N (ν, ω 2 ) for all parameter θ = (ν, ω 2 ) in a compact set Θ ⊆ R × (0, ∞), with true value θ 0 = (ν 0 , ω 2 0 ) ∈ Θ ( Θ is the interior of Θ).

Assumption [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF] The Fisher matrix I(θ 0 ) is nonsingular.

If f (x, ϕ) ∈ C k (R 2 , R), k ≥ 1
, we denote by f (i) the ith derivative with respect to x.

The exact likelihood function

Consider N real valued stochastic processes (X i (t), t ≥ 0), i = 1, . . . , N , with dynamics ruled by [START_REF] Beal | Estimating population kinetics[END_REF]. The processes (W i ) 1≤i≤N and the random variables (r.v.'s) φ i , i = 1, . . . , N are defined on a common probability space (Ω, F , P). Consider the filtration (F t , t ≥ 0) defined by

F t = σ (φ i , W i (s), s ≤ t, i = 1, . . . , N ). As F t = σ (W i (s), s ≤ t) ∧ F i t , with F i t = σ φ i , φ j , W j (s), s ≤ t, j i independent of W i each process W i is a (F t , t ≥ 0)-Brownian motion.
Moreover, the r.v.'s φ i are F 0 -measurable. Under the assumption (H 1 ), the process X i (t) is well defined φ i , X i (t) is strong Markovian adapted to the filtration (F t , t ≥ 0). The N processes φ i , X i (•) 1≤i≤N are independent (this is inherited from the independence of φ i and W i ). For all ϕ ∈ R, x i ∈ R, the following SDE

dX i ϕ (t) = b X i ϕ (t), ϕ + σ X i ϕ (t), ϕ dW i ϕ (t), X i ϕ (0) = x i (2)
admits a unique strong solution process X i ϕ (t), 0 ≤ t ≤ T i adapted to the filtration (F t , 0 ≤ t ≤ T i ). We construct The exact likelihood in the same way as in Delattre ( see [START_REF] Delattre | Maximum likelihood estimation for stochastic differential equations with random effects[END_REF]

[2, Section 3]) Λ N (θ) = N i=1 λ i (X i , θ) , λ i (X i , θ) = R L T i (X i , ϕ)g(ϕ, θ)dϕ (3)
where L T i (X i , ϕ) is the usual conditional likelihood function of the SDE (2). It is worth stressing that is difficult to obtain an exact form of the likelihood in our situation, for two reasons. The first, is that the transition densities of (2) are not explicit in general. The second, is that the marginal densities λ i (X i , θ) are of integral form which cannot be expressed in closed-form, unless we consider a linear model with a specific distribution as in Delattre et al., 2012[2].

Let the processes (X i (t), t ∈ [0, T i ]), i = 1, . . . , N be discretely observed at equidistant time points t i j = jT i n , j = 0, . . . , n and set

∆ i = T i n , X i = X i (t i j ) = X i j , j = 0, . . . , n with T i ≤ T , i = 1, . . . , N
where T 1 , . . . T N , T are fixed. Using transition densities combined with the Markovian nature of ( 2), which the discreet data inherit, we have

Λ N ,n (θ) = N i=1 R n i=1 p ∆ i (X i j , X i j-1 , ϕ)g(ϕ, θ)dϕ (4)
Since the transition densities of (2) are generally not explicit, we propose to use the Gaussian approximation based on Euler scheme. First, we start by standardizing the diffusion functions of X i , i.e., transforming X i into Y i defined as

Y i = γ(X i , ϕ) = X i du σ (u, ϕ) , for allϕ ∈ R.
where any primitive of 1/σ may be selected, i.e.,the constant of integration is irrelevant, because σ (•,

•) ≥ σ 0 > 0 ( otherwise σ (•, •) ≤ -σ 0 < 0 )
. The function γ is increasing (resp decreasing) and invertible for all ϕ ∈ R. applying Ito's Lemma, for i = 1, . . . , N , Y i has unit diffusion, that is :

dY i (t) = µ Y i (t), ϕ dt + dW i (t) , for allϕ ∈ R where µ (y, ϕ) = b γ -1 (y, ϕ), ϕ σ (γ -1 (y, ϕ), ϕ) - 1 2 σ γ -1 (y, ϕ), ϕ
4 Approximating the transition densities we are interested in approximating the transition densities P ∆ i (X i j , X i j-1 , ϕ), i = 1, . . . , N for two adjacent discrete time observations X i j-1 and X i j . To simplify notations, we focus on approximating the transition densities corresponding to X 1 = X. we set W 1 = W , T 1 = T , and ∆ 1 = ∆ = T n . Let P X (∆, x, w) denote the conditional density of X t+∆ = w given X t = x, induced by the model :

dX t = b(X t , ϕ)dt + σ (X t , ϕ)dW t
We define the functions M i (ϕ), i = 1, . . . , N , G(y, ϕ) and Ψ ϕ (∆, x, w) :

M i (ϕ) = ess sup y∈D Y i µ (y, ϕ) + µ 2 (y, ϕ) (5) Ψ ϕ (∆, x, w) = 1 2π∆σ 2 (w, ϕ) exp - 1 2∆ (z -y -∆µ(y, ϕ)) 2 (6)
with y = γ(x, ϕ) and z = γ(w, ϕ)

G(y, ϕ) = y y 0 µ(v, ϕ)dv , for some fixed y 0 ∈ R. (7)
Lemma 1. Under the assumptions ( 1)-( 3), the function M 1 (ϕ) is bounded and P X (∆, x, w) ∼ Ψ ϕ (∆, x, w) as ∆ -→ 0 uniformly in x, w and ϕ.

Proof of Lemma 1 : Let y = γ(x, ϕ), z = γ(w, ϕ) and Y = γ(X, ϕ). we have

µ (y, ϕ) + µ 2 (y, ϕ) ≤ b (x, ϕ) + b(x, ϕ) σ (x, ϕ) σ (x, ϕ) + 1 2 σ (x, ϕ)σ (x, ϕ) + b(x, ϕ) σ (x, ϕ) - 1 2 σ (x, ϕ) 2 ≤ b (x, ϕ) + L 1 2 + σ (x, ϕ) + 2 b 2 (x, ϕ) σ 2 (x, ϕ) + 1 2 σ 2 (x, ϕ) ≤ 2L 2 + L 2 + sup x,ϕ b (x, ϕ) + L σ (x, ϕ) + 1 2 σ 2 (x, ϕ) < ∞ Hence, the functions M i (ϕ), i = 1, . . . ,
P Y (∆, y, z) ∼ e G(z,ϕ)-G(y,ϕ) P W (∆, y, z)
as ∆ -→ 0 uniformly in y, z.

P W (∆, y, z) = d dz P (W t+∆ ≤ z|W t = y) = d dz P(W t+∆ -W t ≤ z -y) = e -1 2∆ (z-y) 2 √ 2π∆
Applying the jacobian formula for the change of density, we get :

P (∆, x, w) = |σ (w, ϕ)| -1 P Y (∆, y, z) ∼ e G(z,ϕ)-G(y,ϕ) 2π∆σ 2 (w, ϕ) × e -1 2∆ (z-y) 2
Hence,

P X (∆, x, w) Ψ ϕ (∆, x, w) ∼ exp G(z, ϕ) -G(y, ϕ) - 1 2∆ (z -y) 2 + 1 2∆ (z -y -∆µ(y, ϕ)) 2 =Ψ ( 8
)
Taking the logarithm of the right hand side in [START_REF] Kessler | Estimation of an ergodic diffusion from discrete observations[END_REF] followed by the absolute value, we have

|log Ψ | ≤ z y µ(u, ϕ)du -µ(y, ϕ)(z -y) + ∆ 2 µ 2 (y, ϕ) ≤ z y [µ(u, ϕ) -µ(y, ϕ)] du + L 2 1 2 ∆ ≤ L 2 1 2 ∆ + sgn(z -y) z y µ(u, ϕ) -µ(y, ϕ) du ≤ L 2 1 2 ∆ + L 2 sgn(z -y) z y u -y du ≤ L 2 1 2 ∆ + L 2 2 (z -y) 2 ≤ L 2 1 2 ∆ + L 2 2 γ(X t+∆ , ϕ) -γ(X t , ϕ) 2 L 2 1 2 ∆ + L 2 2σ 2 0 |X t+∆ -X t | 2 with L 1 := L + 1 2 sup
x,ϕ σ (x, ϕ) and L 2 := sup y,ϕ µ (y, ϕ) . Taking the expectation and using theorem 4 of Gikhman and Skorokhod (1972, p. 48), there exists C > 0 such that,

E |log Ψ | ≤ L 2 1 2 ∆ + CL 2 2σ 2 0 (E(X 2 0 ) + 1) ∆ -→ 0 as ∆ -→ 0
and this holds uniformly in x, w and ϕ.

5 Approximating the likelihood function Λ N ,n (θ)

Remark 2. One may use the approximations given in section 4 instead of the transition densities P ∆ i (X i j , X i j-1 , ϕ). That is not enough to construct a consistent approximate likelihood, since it may not be possible to bound the total asymptotic error of the approximations given by ( 6). Thus we trim the approximations Ψ ϕ (∆ i , X i j-1 , X i j ) by the following compensatory error functions

Φ n,i = exp 1 2n T i 0 µ (Y i (s), ϕ) + µ 2 (Y i (s), ϕ) ds , with Y i (s) = γ(X i (s), ϕ)
We define the approximate likelihood function of the parameter θ by

Λ N ,n (θ) = N i=1 R n j=1 Φ n,i Ψ ϕ (∆ i , X i j-1 , X i j )g(ϕ, θ)dϕ = N i=1 R exp 1 2 U i (ϕ) -1 2 S i (ϕ) (2π∆ i ) n/2 R 1/2 i (ϕ) g(ϕ, θ)dϕ (9)
where

R i (ϕ) = n j=1 σ 2 (X i j , ϕ), U i (ϕ) = T i 0 µ (Y i (s), ϕ) + µ 2 (Y i (s), ϕ) ds and S i (ϕ) = 1 ∆ i n j=1 γ(X i j , ϕ) -γ(X i j-1 , ϕ) -∆ i µ(γ(X i j-1 , ϕ), ϕ) 2 Lemma 3.
Under the assumptions ( 1)-( 3), the approximate likelihood given in ( 9) is consistent as n -→ ∞.

Proof of Lemma 3 :

Let Ψ i n (ϕ) = exp 2 U i (ϕ) -1 2 S i (ϕ) (2π∆ i ) n/2 R 1/2 i (ϕ) for i = 1, . . . , N . We shall show that Ψ i n (ϕ) ∼ n j=1 P ∆ i (X i j , X i j-1 , ϕ) as n -→ ∞. We set y i = Y i (0) = γ(X i (0), ϕ), Y i t = γ(X i (t), ϕ) and F(y) = y y i µ(u, ϕ)du. applying Itô formula to F(Y i t ), we get Y i T i y i µ(u, ϕ)du = T i 0 µ 2 (Y i s , ϕ)ds + 1 2 T i 0 µ (Y i s , ϕ)ds + T i 0 µ(Y i s , ϕ)dW i (s) (10)
For sufficiently large n, (8) yields

n j=1 P ∆ i (X i j , X i j-1 , ϕ) Ψ i n (ϕ) ∼ e F(Y i T i )-1 2 U i (ϕ)-n j=1 µ(Y i j-1 ,ϕ)(Y i j -Y i j-1 )+ 1 2 n j=1 µ 2 (Y i j-1 ,ϕ)∆ i (11)
Taking the logarithm of right hand side in [START_REF] Overgaard | Nonlinear mixed-effects models with stochastic differential equations: implementation of an estimation algorithm[END_REF] and using [START_REF] Liptser | Statistics of random processes I: general theory[END_REF] we get the following term

Σ = 1 2 T i 0 µ 2 (Y i s , ϕ)ds + T i 0 µ(Y i s , ϕ)dW i s - n j=1 µ(Y i j-1 , ϕ)(Y i j -Y i j-1 ) + 1 2 n j=1 µ 2 (Y i j-1 , ϕ)∆ i = - 1 2 T i 0 µ 2 (Y i s , ϕ)ds + T i 0 µ(Y i s , ϕ)dY i s - n j=1 µ(Y i j-1 , ϕ)(Y i j -Y i j-1 ) + 1 2 n j=1 µ 2 (Y i j-1 , ϕ)∆ i = I 1 + 1 2 I 2 with I 1 = T i 0 µ(Y i s , ϕ)dY i s - n j=1 µ(Y i j-1 , ϕ)(Y i j -Y i j-1 ) and I 2 = n j=1 µ 2 (Y i j-1 , ϕ)∆ i - T i 0 µ 2 (Y i s , ϕ)ds
Clearly E(Σ) -→ 0 as n -→ ∞ (see e.g. Bishwal 2008), and the proof is complete.

Consistency of the AMLE

Consider x i = x and T i = T , for i = 1, . . . , N so that the set-up boils down to the iid situation. Following Schervish 1995[15, Theorem 7.49 and Lemma 7.54], we propose to prove sufficient conditions ensuring consistency of the Maximum Likelihood Estimator (MLE) of the parameter θ based on Λ N ,n . Let C T 1 denote the space of real continuous functions (x(s) : s ∈ [0, T 1 ]), defined on [0, T 1 ], endowed with σ -field B 1 associated with the topology of uniform convergence [0, T 1 ]. Consider the sequence X i ∞ i=1 of random variables conditionally iid given θ, with density function λ 1 (x, θ) on the space (C T 1 , B 1 ).

λ 1 (x, θ) ∝ R exp 1 2 U 1 (ϕ) -1 2 S 1 (ϕ) R 1 (ϕ) 1/2 g(ϕ, θ)dϕ. ( 12
)
Fix θ 0 ∈ R × (0, ∞), and define for each

M ⊆ R × (0, ∞) and x ∈ C T 1 Z (M, x) = inf θ∈M log λ 1 (x, θ 0 ) λ 1 (x, θ)
Lemma 4. Under the assumptions ( 1)-( 4), and for each θ θ 0 , there exists an open set N θ such that θ ∈ N θ and

E θ 0 Z(M, X i ) > -∞.
Proof of Lemma 4 : Clearly U 1 (ϕ) is bounded. Indeed, for all ϕ ∈ R, we have

U 1 (ϕ) ≤ (L 2 1 + L 2 )T 1 2S 1 < ∞. Let g(ϕ, θ)dϕ = N (ν, ω 2 ), θ = ω 2 ). For θ θ 0 set A = log R exp 1 2 U 1 (ϕ) -1 2 S 1 (ϕ) R 1 (ϕ) 1/2 e -1 ω 2 0 (ϕ-ν 0 ) 2 dϕ and B = log R exp 1 2 U 1 (ϕ) -1 2 S 1 (ϕ) R 1 (ϕ) 1/2 e -1 ω 2 (ϕ-ν) 2
dϕ. To prove lemma 4 it suffices to bound A and -B from below by random variables with finite expectation. Applying the fact that for all a, b ∈ R; (a + b) 2 ≤ 2(a 2 + b 2 ), we can assert that

S 1 (ϕ) = 1 ∆ 1 n j=1 Y 1 j -Y 1 j-1 -∆ 1 µ(Y 1 j-1 , ϕ) 2 ≤ 2 ∆ 1 n j=1 (Y 1 j -Y 1 j-1 ) 2 + 2∆ 1 n j=1 µ 2 (Y 1 j-1 , ϕ) ≤ 2 ∆ 1 σ 2 0 n j=1 (X 1 j -X 1 j-1 ) 2 S 2 + 2T 1 L 2 1 S 3 Therefore A = log R exp 1 2 U 1 (ϕ) -1 2 S 1 (ϕ) R 1 (ϕ) 1/2 e -1 ω 2 0 (ϕ-ν 0 ) 2 dϕ ≥ -(S 1 + S 2 + S3) + log R e -1 ω 2 0 (ϕ-ν 0 ) 2 R 1 (ϕ) 1/2 dϕ
Let K be a linear growth constant for σ , i.e., σ 2 (x, ϕ)

≤ K 2 (1 + x 2 + ϕ 2 ) , for all x, ϕ ∈ R R 1 (ϕ) = n j=1 σ 2 (X 1 j , ϕ) ≤ K 2n 1 + sup 0≤j≤n (X 1 j ) 2 + ϕ 2 n ≤ K 2n 2 n-1 2 + sup 0≤j≤n (X 1 j ) 4n + ϕ 4n =⇒ R 1/2 1 (ϕ) ≤ K n 2 (n-1)/2 √ 2 + sup 0≤j≤n (X 1 j ) 2n + ϕ 2n ≤ K n 2 (n-1)/2 √ 2 + sup 0≤j≤n (X 1 j ) 2n 1 + ϕ 2n ≤ K n 2 (n-1)/2 √ 2 + sup 0≤j≤n (X 1 j ) 2n e 2nϕ Hence, R e -1 ω 2 0 (ϕ-ν 0 ) 2 R 1 (ϕ) 1/2 dϕ ≥ K -n 2 (1-n)/2 √ 2 + sup 0≤j≤n (X 1 j ) 2n -1 R e -1 ω 2 0 (ϕ-ν 0 ) 2 -2nϕ dϕ Taking the logarithm yields log R e -1 ω 2 0 (ϕ-ν 0 ) 2 R 1 (ϕ) 1/2 dϕ ≥ -log(K n 2 (n-1)/2 ) -log √ 2 + sup 0≤j≤n (X 1 j ) 2n + log R e -1 ω 2 0 (ϕ-ν 0 ) 2 -2nϕ dϕ ≥ -log(K n 2 (n-1)/2 ) + log R e -1 ω 2 0 (ϕ-ν 0 ) 2 -2nϕ dϕ - √ 2 -sup 0≤j≤n (X 1 j ) 2n
By virtue of Delattre et al., 2012 [2, Proposition 1] and the fact that all moments of φ i are finite, the random variables S 2 and sup 0≤j≤n (X 1 j ) 2n have finite expectation

under P θ 0 , thereby E θ 0 (A) > -∞. It remains to show that E θ 0 (-B) > -∞. Note that e 1 2 U 1 (ϕ) ≤ e (L 2 1 +L 2 )T 1 /2 = e S 1 and 1
R 1/2 1 (ϕ) ≤ σ -n 0 .
Simple computations leads to the following inequality

-B ≥ S 1 + log σ -n 0 R e -1 ω 2 (ϕ-ν) 2 dϕ
From the computations given above, it follows that E θ 0 Z(M, X i ) > -∞ and N θ could be any neighborhood of θ. We claim that the function θ → λ 1 (x, θ) is continuous. In fact the smoothness of g(ϕ, •) on the compact set

Θ = [µ, µ] × [ω 2 , ω 2 ]
for any fixed ϕ ∈ R provides the lipschitz criterion, thus we construct a density function h(.), so that

|g(ϕ, θ 1 ) -g(ϕ, θ 2 )| ≤ C Θ h(ϕ) θ 1 -θ 2 , ∀(θ 1 , θ 2 ) ∈ Θ,
and λ 1 (x, θ 1 )λ 1 (x, θ 2 ) ≤ MC Θ σ -n 0 e S 1 θ 1θ 2 where C Θ is some positive constant depending only on the parameter space Θ and M is the positive proportionality constant in [START_REF] Picchini | Stochastic differential mixed-effects models[END_REF], hence the continuity of λ 1 (x, θ) at θ for every θ, a.s P θ 0 . Then all conditions of Lemma 7.54 in Schervish 1995 [START_REF] Schervish | Theory of Statistics[END_REF] are fulfilled. Hence the following statement is fully established.

Theorem 5. The AMLE θ N ,n of θ is almost surely consistent under P θ 0 , that is

θ N ,n P θ 0 a.s -→ N →∞ θ 0 .
7 Asymptotic normality of the approximate maximum likelihood estimator

Let l N ,n (θ) be the logarithm of the approximate likelihood defined as :

l N ,n (θ) = log Λ N ,n (θ) ∝ N i=1 R exp 1 2 U i (ϕ) -1 2 S i (ϕ) R i (ϕ) 1/2 g(ϕ, θ)dϕ
The smoothness of the function (ϕ, θ) → g(ϕ, θ) allowed us to compute the score function and the Fisher Information matrix. Here we give the exact formula of the score function with integral components with finite expectation under P θ 0 . Let θ = (θ 1 , θ 2 ) = (ν, ω 2 ). Note that there exist functions P sk (ϕ, θ), P s (ϕ, θ), k, s = 1, 2 polynomial in ϕ such that ∂ θ s g(ϕ, θ) = P s (ϕ, θ)g(ϕ, θ), ∂ 2 θ s g(ϕ, θ) = P ss (ϕ, θ)g(ϕ, θ) and ∂ θ s ∂ θ k g(ϕ, θ) = P sk (ϕ, θ)g(ϕ, θ). For every function f (ϕ, θ), we define the functions f i (θ) as

f i (θ) = R exp 1 2 U i (ϕ) -1 2 S i (ϕ) R i (ϕ) 1/2 f (ϕ, θ)g(ϕ, θ)dϕ
The functions f i (θ) are well defined if f (ϕ, θ) is polynomial in ϕ. Indeed, for any θ

E f i (θ) ≤ E R exp 1 2 U i (ϕ) -1 2 S i (ϕ) R i (ϕ) 1/2 |f (ϕ, θ)| g(ϕ, θ)dϕ ≤ e (L 2 1 +L 2 )T i /2 σ n 0 R |f (ϕ, θ)| g(ϕ, θ)dϕ < ∞ Let P (ϕ, θ) = 1.
Simple computations lead to the following form of the score function

S(ν, σ 2 ) = N i=1 P i 1 (θ) P i (θ) -1
N i=1 P i 2 (θ) P i (θ) -1 and the Fisher Information matrix given by:

I(ν 0 , σ 2 0 ) = - E θ 0 Σ 1.1 E θ 0 Σ 1.2 E θ 0 Σ 2.1 E θ 0 Σ 2.2 ,
where

Σ 1.1 = P 1 11 (θ 0 ) -P 1 1 (θ 0 ) 2 P 1 (θ 0 ) -2 , Σ 2.2 = P 1 22 (θ 0 ) -P 1 2 (θ 0 ) 2 P 1 (θ 0 ) -2
and Σ 1.2 = P 1 12 (θ 0 ) -P 1 1 (θ 0 ) P 1 2 (θ 0 ) P 1 (θ 0 ) -2 . To prove the finiteness of the Fisher matrix, it remains to show that P 1 (θ 0 ) -2 is bounded from above by a random variable with finite expectation. For that purpose we add the following assumption Assumption (6) The processes X i , i = 1, . . . , N are discretely observed at n 2 + 1 equidistant time points t i j = jT i /n 2 , j = 0, . . . , n 2 ,

∆ i = T i /n 2
It is well known that under the assumption (6),

n 2 j=1 W 1 (t j ) -W 1 (t j-1 ) 2 -→ T 1
almost surely under P θ 0 . Fix ε > 0, there exist n 0 ∈ N, such that

n 2 j=1 W 1 (t j ) -W 1 (t j-1 ) 2 ≤ T 1 + ε for all n ≥ n 0 Hence, P 1 (θ 0 ) ≥ e -S 1 -(T 1 +ε)/2∆ 1 K n 2 ( √ 2 + sup j (X 1 j ) 2n 2 ) R e -2n 2 ϕ g(ϕ, θ)dϕ C =⇒ P -1 (θ 0 ) ≤ e S 1 +(T 1 +ε)/2∆ 1 C K n 2 ( √ 2 + sup j (X 1 j ) 2n 2 ) =⇒ P -2 (θ 0 ) ≤ C 1 + sup j (X 1 j ) 4n 2
The right hand side of the last inequality is of finite expectation. Hence the Fisher matrix is well defined. Proposition 6. Under the assumptions ( 1)-( 6), the MLE based on Λ N ,n 2 (θ) is asymptotically normal.

Proof of Proposition 6 : The proof is based on the result given in Schervish 1995 [15, theorem 7.63]. We state first this result and check their regularity conditions.

Lemma 7. [START_REF] Schervish | Theory of Statistics[END_REF]) Let Θ be a subset of R p , and let {X N } ∞ N =1 be conditionally iid given θ each with density f X 1 |θ (.|θ). Let θk be an MLE. Assume that θ N P -→ θ under P θ for all θ. Assume that f X 1 |θ (x|θ) has continuous second partial derivatives with respect to θ and that differentiation can be passed under the integral sign. Assume that there exists H r (x, θ) such that, for each θ 0 ∈ Int(Ω) and for each k, j,

sup θ-θ 0 ≤r ∂ 2 ∂ θ k ∂ θ j log f X 1 |θ (x|θ 0 ) - ∂ 2 ∂ θ k ∂ θ j log f X 1 |θ (x|θ) ≤ H r (x, θ 0 ); ( 13
)
with lim r→0 E θ 0 H r (X, θ 0 ) = 0. ( 14)

Assume that the Fisher Information matrix is finite and nonsingular. Then, under

P θ 0 , √ N θ N -θ 0 D =⇒ N 0, I -1 X 1 (θ 0 ) . ( 15
)
In section 6, we proved almost sure consistency of the ALME θ N ,n . Hence θ N ,n 2 P -→ θ for all θ. We have shown that the logarithm of the approximate likelihood l N ,n (θ) has second partial derivatives with finite expectation (section 7). We can also prove, in similar fashion that l N ,n (θ) has third partial derivatives with finite expectation. Hence ( 13) and ( 14) clearly hold. The Fisher matrix I(θ 0 ) is nonsingular by hypothesis [START_REF] Donnet | Parametric inference for mixed models defined by stochastic differential equations[END_REF]. Therefore, asymptotic normality of the MLE based on Λ N ,n 2 (θ), of the form [START_REF] Schervish | Theory of Statistics[END_REF] holds in this case.

  N are bounded. By corollary 2 in Downes (2008), we have


S Shimasaki
email: shimasaki@t.kagawa-nct.ac.jp

S Taniguchi

Separation Efficiency of Inclusion Particles from Liquid Aluminum by Electromagnetic Cyclone Separator

Keywords: electromagnetic cyclone separator, inclusion particles, liquid cyclone, electromagnetic separation

A new inclusion separation device has been proposed; it is the electromagnetic cyclone separator, which is composed of a liquid cyclone and an electromagnetic induction coil. Fluid flow and separation efficiency of the electromagnetic cyclone separator were investigated by cold model experiments. CFD analysis was conducted to predict flows and pressure drop of the separator and compared with the experimental results. A mathematical model based on analysis of particle trajectory was formed to estimate separation efficiency of the device and predicted performance agreed with the observation.

Introduction

It is well known that inclusion particles, which remain in final products, are very harmful for material properties. In order to remove inclusion particles from molten metal, numerous separation ways have been proposed for industrial refining processes. Turchin et al. firstly applied a cyclone concept to impurities removal from molten aluminum [1]. The concept has been followed by Courtenay's research group [2] due to attractive characteristics of the device such as a high-throughput of molten aluminum, the absence of moving parts and low operational costs etc.

We have proposed a new separation device for removal of inclusion particles from molten aluminum; it is the electromagnetic cyclone separator, which is composed of a liquid cyclone and an electromagnetic induction coil [START_REF] Nakabayashi | 7th Int. Conf. on Electromagnetic Processing of Materials (EPM2012)[END_REF]. Fig. 1 shows the schematic diagram of the electromagnetic cyclone separator. The liquid cyclone is a device to separate inclusion particles from molten aluminum on the basis of a density difference between the particle and the suspension medium. Inclusion particles are also separated by an electromagnetic separation force [START_REF] Shimasaki | Proc. of 13th Int. Conf. on Aluminum Alloys (ICAA13)[END_REF], which is derived from a difference of electric conductivity between them. It is expected that the two separation forces will separate inclusion particles complementary to each other. In this study, cold model experiment was conducted in order to investigate characteristics of fluid flow of the cyclone. Experimental results such as velocity profiles and pressure drops were compared with previous studies and CFD analysis. Separation efficiency of the electromagnetic cyclone separator was estimated from the mathematical model, which is based on analysis of particle trajectory. The predicted performance of the cyclone was compared with experimental results of cold mode experiments, in which the magnetic separation was used instead of the electromagnetic separation.

Experiments

Schematic diagram of the cold model experiment is shown in Fig. 2. Filtered tap water was pumped up to the upper tank and a level of the tank was maintained constantly by the overflow. The valve and the flow meter were used to regulate a flow rate of inlet water. Small amount of suspension was injected by the syringe for separation tests. Temperature of the circulating water was kept constant T = 1 ºC by blocks of ice floating in the tanks. Two cyclones of different size were used to evaluate size effects.

Magnetic particles and the sheeted rubber magnet were used in the cold model experiment in order to mimic the electromagnetic separation. Electromagnetic separation cannot be applied to suspension with liquid medium of tap water, because the separation is valid only for electrical conductive fluids. The magnetic particles are made of ferrite (spherical shape, density: 4900 kg•m -3 , diameter: 25-65 μm) and magnetic force acting on the particles was obtained by preliminary experiments as a function of distance from the magnet sheet.

Characteristics of fluid flow in the cyclone were observed by the following methods: flow visualization with coloured ink, analysis of residence time distribution, velocity field taken by PIV, measurement of pressure drop by the cyclone. Separation efficiencies of the cyclone were also measured with/without the magnet sheet for various condition of flow rate, cyclone size and particle diameter.

Flow characteristics of the liquid cyclone

Fig. 3 shows tangential velocity of the cyclone obtained by PIV measurement. The result shows a typical flow profile of liquid cyclone [START_REF] Bergstöm | [END_REF]; the tangential velocity v θ increased from the sidewall toward the center, reached a maximum value and then suddenly decreased. The velocity seems to be proportional to the flow rate. A centrifugal force acting on a particle in the liquid cyclone is proportional to v θ 2 /r, the separation force in the center region is much higher than that in the vicinity of the side wall. Results of CFD calculation are also plotted on Fig. 2. Axially symmetric calculations were carried out in order to reduce calculation costs; the three components of flow (v r , v θ , v z ) were calculated with the two coordinates (r, z). The numerical evaluation gave somewhat higher values than the observation (10-20 % higher), however, the overall trend was predicted well. Pressure drop of the liquid cyclone is one of the most important operation conditions. Fig. 4 shows the relation between hydraulic head loss and flow rate for two different size of the cyclone. The hydraulic head loss increased with increasing of the flow rate. The smaller cyclone had the higher pressure drop. Pressure drop could be estimated by the empirical equation proposed by Dahlstrom [6]:

Q h = c D in ⋅ D out ( ) 0.9 (1)
where Q is a flow rate, c is a constant depending on geometry of the cyclone, D in and D out are diameters of the inlet and the outlet of the cyclone respectively. It can be seen from the figure that the agreements between the experimental results and Dahlstrom's equation were quit well. Pressure drops evaluated by CFD calculation were also agrees well.

Separation efficiency

Fig. 5 shows the relation between separation efficiency and particle diameter for various conditions. A mathematical model was used to predict the removal efficiency. It was based on an analysis of particle trajectory, and the detail of the model was explained on the literatures [START_REF] Nakabayashi | 7th Int. Conf. on Electromagnetic Processing of Materials (EPM2012)[END_REF]7]. Larger particle was easer to remove, because both of the separation forces (centrifugal and magnetic forces) were proportional to the volume of the particle. We can see that the removal efficiency was improved by applying the magnetic separation, especially for the case of 15 L/min. As mentioned above, the centrifugal force highly depends on the tangential velocity, the relative importance of each separation force varied with the flow rate. The mathematical model succeeds to predict the removal efficiency semi-quantitatively. It is expected that the model can be used to predict a performance of the electromagnetic cyclone separator in aluminum refining process.

Conclusions

Flow characteristics of liquid cyclone such as velocity profile and pressure drop were investigated by the experiments and the numerical evaluations. Separation efficiencies with/without magnetic separation were also investigated. It was found that the removal efficiency of the cyclone was effectively improved by applying magnetic separation.

Fig. 1 :

 1 Fig. 1: Schematic diagram of the electromagnetic cyclone separator.

Fig. 2 :

 2 Fig. 2: Schematic diagram of the electromagnetic cyclone separator.

1 Fig. 3 :

 13 Fig. 3: Profiles of tangential velocity of the liquid cyclone.

Fig. 4 :

 4 Fig. 4: Relation between head loss and flow rate.

5 :

 5 only Exp. cyclone + magnet Cal. cyclone only Cal. cyclone + magnet particle diameter, d p / μm removal efficiency, η / % only Exp. cyclone + magnet Cal. cyclone only Cal. cyclone + magnet particle diameter, d p / μm removal efficiency, η / % (a) 15 L/min (b) 30 L/min Fig. Relation between head loss and flow rate.

Acknowledgment

This work was financially supported by JSPS KAKENHI, Grant-in-Aid for Scientific Research (A) Grant Number 22246097.