Malik Saad

Mohammed Missen

Mickaël Attik

Antoine Coustaty

Cyril Doucet

Faucher

SentiML++ : An Extension of the SentiML Sentiment Annotation Scheme

In this paper, we propose SentiML++, an extension of SentiML with a focus on annotating opinions answering aspects of the general question "who has what opinion about whom in which context?". A detailed comparison with SentiML and other existing annotation schemes is also presented. The data collection annotated with SentiML has also been annotated with SentiML++ and is available for download for research purpose.

Introduction

The semantic annotation of opinions is one of the very important tasks of opinion mining. Semantic annotations are very important both for training machine learning approaches and for evaluating opinion mining methods. Unfortunately, there have been hardly any serious proposal attempts of appropriate annotation schemas until recently when SentiML [START_REF] Di Bari | Sentiml: Functional annotation for multilingual sentiment analysis[END_REF], OpinionMiningML [START_REF] Robaldo | Opinionmining-ml[END_REF] and EmotionML [START_REF] Schröder | Emotionml -an upcoming standard for representing emotions and related states[END_REF] were proposed. In this paper, we discuss, compare and identify the positives and negatives of these annotation schemes. Following this overview, we propose Sen-tiML++, an extension of SentiML that addresses several shortcomings of the state of the art.

SentiML The SentiML annotation schema [START_REF] Di Bari | Sentiml: Functional annotation for multilingual sentiment analysis[END_REF] follows a conventional sentiment annotation style and is based on Appraisal Framework (AF) [START_REF] Martin | The language of evaluation, Appraisal in English[END_REF] which is a strong linguistically-grounded theory. AF helps to define appraisal types (affect, judgments and appreciation) within the modifier tag which is another positive point to be noted in SentiML. With a very simple annotation scheme, SentiML is popular because adopting its annotation scheme does not require to acquire any specific skills. However, concerns can be raised about SentiML.

OpinionMiningML OpinionMiningML [START_REF] Robaldo | Opinionmining-ml[END_REF] is an XML-based formalism that allows tagging of attitude expressions for features or objects as found in a textual segment. It targets extraction of feature-based opinion expressions but its scope is limited to proposing an annotation schema. Besides this, the structure of OpinionMiningML is not straightforward and can be threatened by challenges for feature and relation extraction while developing an automatic tagger for this annotation scheme.

EmotionML EmotionML [START_REF] Schröder | Emotionml -an upcoming standard for representing emotions and related states[END_REF] aims to make concepts from major emotion theories available in a broad range of technological contexts. Being informed by the effective sciences, EmotionML recognises the fact that there is no single agreed representation of effective states, nor of vocabularies to use. Therefore, an emotional state < emotion > can be characterised using four types of descriptions: < category >, < dimension >, < appraisal > and < action -tendency >. Furthermore, the vocabulary used can be identified.

SentiML example Throughout the article, we will use the following sentence as a running example: "The U.S. State Department on Tuesday (KST) rated the human rights situation in North Korea "poor" in its annual human rights report, casting dark clouds on the already tense relationship between Pyongyang and Washington." Relevant annotations in SentiM L are given below:

<APPRAISALGROUP id="A0" fromID="T0" fromText="situation" toID="M0" toText="poor" orientation="negative "/> <APPRAISALGROUP id="A1" fromID="M1" fromText="dark" toID="T1" toText="clouds" orientation="negative"/ > <APPRAISALGROUP id="A2" fromID="M2" fromText="tense" toID="T2" toText="relationship" orientation=" negative"/> <MODIFIER id="M0" start="201" end="205" text="poor" attitude="appreciation" orientation="negative" force=" normal" polarity="unmarked"/> <MODIFIER id="M1" start="250" end="254" text="dark" attitude="appreciation" orientation="negative" force=" normal" polarity="unmarked"/> <MODIFIER id="M2" start="277" end="282" text="tense" attitude="appreciation" orientation="negative" force=" normal" polarity="unmarked"/> <TARGET id="T0" start="175" end="184" text="situation" type="thing" orientation="neutral"/> <TARGET id="T1" start="255" end="261" text="clouds" type="thing" orientation="ambiguous"/> <TARGET id="T2" start="283" end="295" text="relationship" type="thing" orientation="neutral"/> OpinionMiningML example below annotated using the OpinionM iningM L syntax:

<COMMENT id="1" ontologyreference="1"> <FRAGMENT id="1"> The U.S. State Department on Tuesday (KST) rated the human rights situation in North Korea "poor" </FRAGMENT> <FRAGMENT id="2"> in its annual human rights report, casting dark clouds on the already tense relationship between Pyongyang and Washington.</FRAGMENT> <APPRAISAL polarity="negative" intensity="medium"> <FACETREFERENCE>1</FACETREFERENCE> <FRAGMENTREFERENCE>1</FRAGMENTREFERENCE> <FRAGMENTREFERENCE>2</FRAGMENTREFERENCE> </APPRAISAL> <APPRAISAL polarity="negative" intensity="medium"> <FACETREFERENCE>2</FACETREFERENCE> <FRAGMENTREFERENCE>1</FRAGMENTREFERENCE> <FRAGMENTREFERENCE>2</FRAGMENTREFERENCE> </APPRAISAL> </COMMENT> EmotionML example presented hereby an example annotated using the EmotionML syntax: <emotionml xmlns="http://www.w3.org/2009/10/emotionml" xmlns:meta="http://www.example.com/metadata" category-set="http://www.w3.org/TR/emotion-voc/xml#occ-categories"> <info> <meta:doc>Example taken from annotation of SentiML </meta:doc> </info> The U.S. State Department on Tuesday (KST) <emotion> < category name="reproach"/> rated the human rights situation in North Korea "poor" </emotion> in its annual human rights report <emotion> <category name="disappointment"/> casting dark clouds </emotion> on the already <emotion> < category name="disappointment"/> tense relationship between Pyongyang and Washington. </emotion> </emotionML>

Comparison

In this section we give a comparison of annotations schemes from different perspectives, summarized in Table 1. From this comparison, it can be concluded that SentiML has a larger scope and it is equipped with a more affordable vocabulary with respect to the previous work [START_REF] Cambria | New avenues in opinion mining and sentiment analysis[END_REF][START_REF] Liu | A survey of opinion mining and sentiment analysis[END_REF][START_REF] Missen | Opinion mining: reviewed from word to document level[END_REF][START_REF] Pang | Opinion mining and sentiment analysis[END_REF]. Hence, we find it the most suitable choice for our current work.

SentiML ++

In this section, we provide an extension of SentiML considering the work of Bing Liu [START_REF] Liu | Sentiment analysis and subjectivity[END_REF] as a reference and find out that most of the aspects defining an opinion seem to be missing (completely or partially). For example, SentiML works on sub-sentence level and hence leaves actual holder and target entities of sentiment of a sentence unmarked. As far as opinion orientation is concerned, it deals with prior polarities in a better way than the contextual ambiguities. It does not recognize the topic of a sentence, hence fails to identify topic-based contextual ambiguities. Similarly, the contexts defined by cultural phrases and emoticons cannot be identified using SentiML. Flexibility and completeness are important characteristics of an annotation scheme [START_REF] Oren | What are semantic annotations? Tech. rep[END_REF] Scope EmotionML, the W3C standard, is an effort to cover all aspects of emotions globally in all concerned fields while OpinionMiningML and SentiML limit themselves to domains of IR and NLP. In our view, SentiML has larger scope as compared to OpinionMiningML which limits itself only to feature-based sentiment analysis. Complexity Because of its larger scope than the other two annotation schemes, EmotionML can be considered complex and less user-friendly. SentiML is much easier to use than OpinionMiningML because of the vocabulary of its annotation scheme which matches the vocabulary used in research work of this field, (i.e. concepts like holder, target, modifier, etc.). Vocabulary EmotionML text annotation is equipped with new and broader vocabulary while targets of SentiML and OpinionMiningML are more specific. SentiML revolves around the concept of modifier and targets of the sentiments while OpinionMiningML's focus remains on sentiment relevant to features of the objects and equips itself with numerous meta-tags.

Structure

The structure of all three annotation schemes follow XML based format but OpinionMiningML defines granularity better than the other annotation schemes by reaching up to the feature level. Contextual Ambiguities While all three annotation schemes focus on defining semantics of the expression types, such as appreciation, suggestion etc., SentiML also tackles the issue of contextual ambiguities which is a big research challenge in the field of opinion mining. Theoretical Grounds EmotionML is a W3C standard and has been recommended after years of discussion and debate between experts and stakeholders while SentiML is based on the Appraisal Framework, a strong linguistically-grounded theory. OpinionMiningML, in the contrary, has not been reported to have a theoretical basis. Granularity Level Another distinction that can be observed while comparing these annotation schemes is their granularity level. EmotionML is found to be operating on the sentence level while OpinionMiningML and SentiML are rather focusing on sub-sentence levels. While it is mandatory to analyze on a sub-sentence (or word level) to find correct sentiments, the sentence is however a more logical unit of discourse. CompletenessCompleteness [START_REF] Oren | What are semantic annotations? Tech. rep[END_REF] is one of the most important features of annotation schemes. Completeness demands whether an annotation deals with all or most of the real world scenarios. Observing all of these annotations i.e. SentiML, EmotionML and OpinionMiningML in light of this particular characteristic, no annotation scheme seems to satisfy it. SentiML and EmotionML behave similarly by just annotating sentimental expressions while OpinionMiningML goes a step forward by targeting features of objects but fails to deal with contextual aspects of opinions, one of the major problems of opinion mining. Flexibility EmotionML, knowing the challenges of different domains, gives freedom to use whatever vocabulary of emotional states suited to a particular domain. OpinionMiningML follows EmotionML in this regard and gives the choice of building a domain-based ontology as required, a possibility that SentiML lacks.

Table 1. Comparison of Annotation Schemes

and unfortunately, SentiML fails to have both of these characteristics. Identification of opinion words and their polarity with respect to a given topic could help resolving contextual ambiguities. Therefore, we propose to take topic identification into account in SentiML++ by proposing <TOPIC> element. A good share of the opinions generally found on the web are expressed informally. This includes the use of emoticons and sarcastic phrases or cultural expressions (e.g., "bored to death", "dressed to kill", etc.) that could invert the semantics of the text surrounding them. Identification of such contexts could aid the automated detection of such opinion inversions.In SentiML++, we deal with this problem using <INFORMAL> element.

SentiML++ Example SentiML++ operates on both levels i.e. phrase as well as sentence level.

In this section, we annotate the same sentence as an example using SentiML++. The annotation includes the introduction of sentence-level markups like <HOLDER>, <TARGET>, <ORIENTATION>, <TOPIC> and <INFORMAL> while only <HOLDER> is introduced on the sub-sentence level. All of these markups come under the main markup of <SENTENCE> while sub-sentence level annotation comes under <PHRASES>. The markup <SENTENCE> includes one attribute called "type" with possible values "general" or "informal". The value "general" is used for ordinary sentences while the value "informal" is only used when the whole identified sentence is an idiom, a metaphor or an emoticon. When the value "informal" is used, only the <INFORMAL> markup plays its role while other markups are discarded because they are rendered meaningless. The <HOLDER> markup was introduced on the sub-sentence level to identify the holders even at this smaller granularity, if needed.

<SENTENCE ID=0 type="general"> <HOLDER id="H0" start="33" end="53" text="U.S. State Department" type="organization" orientation="neutral "/> <TARGET id="T0" start="108" end="117" text="North Korea" type="country" orientation="neutral"/> <ORIENTATION polarity="negative"> <TOPIC id="1" ref="www.dmoz.org" chain="Society:Issues:Human Rights and Liberties"> <PHRASES> <APPRAISALGROUP id="A0" fromID="T0" fromText="situation" toID="M0" toText="poor" Holder="H0" orientation="negative"/> <APPRAISALGROUP id="A1" fromID="M1" fromText="dark" toID="T1" toText="clouds" Holder="H1" orientation="negative"/> <APPRAISALGROUP id="A2" fromID="M2" fromText="tense" toID="T2" toText="relationship" Holder="H2" orientation="negative"/> <MODIFIER id="M0" start="201" end="205" text="poor" attitude="appreciation" orientation="negative" force= "normal" polarity="unmarked"/> <MODIFIER id="M1" start="250" end="254" text="dark" attitude="appreciation" orientation="negative" force= "normal" polarity="unmarked"/> <MODIFIER id="M2" start="277" end="282" text="tense" attitude="appreciation" orientation="negative" force= "normal" polarity="unmarked"/> <TARGET id="T1" start="175" end="184" text="situation" type="thing" orientation="neutral"/> <TARGET id="T2" start="255" end="261" text="clouds" type="thing" orientation="ambiguous"/> <TARGET id="T3" start="283" end="295" text="relationship" type="thing" orientation="neutral"/> <HOLDER id="H1" start="33" end="53" text="U.S. State Department" type="organization" orientation="neutral "/> <HOLDER id="H2" start="33" end="53" text="U.S. State Department" type="organization" orientation="neutral "/> <HOLDER id="H3" start="33" end="53" text="U.S. State Department" type="organization" orientation="neutral "/> </PHRASES> </SENTENCE> It must be noted that <APPRAISALGROUP> in SentiML++ links modifier, target and holder identified at the sub-sentence level. This is in contradiction with SentiML [START_REF] Di Bari | Sentiml: Functional annotation for multilingual sentiment analysis[END_REF], where only modifier and target are linked. <Holder> and <Target> elements are found on both levels i.e. sentence and sub-sentence level. Natural language processing techniques such as syntactic parsing can be helpful in identifying these elements on both granularity levels.

Conclusions and Future Work

In this paper, we proposed SentiML++, an extension of SentiML. We proposed to add target (on sentence level), holder, topic and informal sentence identification as part of SentiML++. SentiML++ adds flexibility to SentiML by giving freedom of choice for a taxonomy when annotating the topic of a sentence. As part of our future work, we plan to further enhance SentiML++ by modeling relations between its elements. The idea is to propose a more suitable model for the semantic web so that state-of-the-art semantic web tools can be leveraged to exploit its semantics.


Nki Echenim

Frederique Clément

Michel Sorine

Multi-scale modeling of follicular ovulation as a reachability problem

Keywords: biomathematics, conservation laws, method of characteristics, control theory, backwards reachable sets

Multi-scale modeling of follicular ovulation as a reachability problem

Nki Echenim, Frédérique Clément, Michel

Introduction

The development of ovarian follicles is a crucial process for reproduction in mammals, as its biological meaning is to free fertilizable oocyte(s) at the time of ovulation. A better understanding of follicular development is both a clinical and zootechnical challenge; it is required to improve the control of anovulatory infertility in women, as well as ovulation rate and ovarian cycle chronology in domestic species. Within all the developing follicles, very few actually reach the ovulatory size; most of them undergo a degeneration process, known as atresia [START_REF] Greenwald | Follicular development and its control[END_REF]. The ovulation rate (number of ovulatory follicles per cycle) results from an FSH (Follicle Stimulating Hormone)-dependent follicle selection process. FSH acts on the cells surrounding the oocyte, the granulosa cells, and controls both their commitment towards either proliferation, differentiation or apoptosis and their ability to secrete hormones such as estradiol. The whole estradiol output from the ovaries is responsible for exerting a negative feedback on FSH release. Following the subsequent fall in plasmatic FSH levels, most of the follicles undergo atresia and only the ovulatory ones survive in the FSH-poor environment. We have proposed a mathematical model, using both multi-scale modeling and control theory concepts, to describe the follicle selection process [START_REF] Echenim | Multi-scale modeling of the follicle selection process in the ovary[END_REF]. For each follicle, the cell population dynamics is ruled by a conservation law with variable coefficients, which describes the changes in age and maturity of the granulosa cell density. A control term, representing FSH signal, intervenes both in the velocity and loss terms of the conservation law. Two acting controls are distinguished: a global control resulting from the ovarian feedback and corresponding to FSH plasmatic levels, and a local control, specific to each follicle, accounting for the modulation in FSH bioavaibility due to follicular vascularization. Both ovulation triggering and follicular ovulation depend on the reaching of a target. In this paper, we aim at using control theory to characterize follicular trajectories and the control laws leading to ovulation or atresia. The macroscopic phenomenon of follicular ovulation is considered as a reachability problem for the microscopic characteristics associated with the conservation law. The problem is enounced in details and solved in open-loop. Follicles are assumed to ovulate if their state variables, (age, maturity and cell density) reach a given target set. Backwards reachable sets theory is used to define the initial conditions compatible with latter ovulation, for a given control panel. In section 2, the conservation laws describing the model for follicle selection and their characteristics are presented. In section 3, the target corresponding to follicular ovulation is defined. Reachable sets theory is used to solve the corresponding control problem and simulation results are discussed. Section 4 is devoted to the discussion and perspectives.

2 Follicle selection model 2.1 Controlled conservation laws for granulosa cells Following [START_REF] Echenim | Multi-scale modeling of the follicle selection process in the ovary[END_REF], cells are characterized by their positions within or outside the cell cycle and by their sensitivities to FSH. This leads to distinguish 3 cellular phases within the granulosa cell population. Phases 1 and 2 correspond to the proliferation phases (describing respectively the G1 phase and S to M phases of the cell cycle), and phase 3 corresponds to the differentiation phase, after cells have exited the cell cycle (see Figure 1). More precisely, the position of a cell at time t ∈ [0, T ] (T > 0) is defined by its age and maturity. The cell age a is a marker of progression within the cell cycle (in phases 1 and 2) and evolves as time t outside the cycle (in phase 3). The duration of phase 1 is a 1 > 0 and the total cycle duration is a 2 (so that phase 2 duration is a 2 -a 1 > 0). The maturity marker γ is used to sort the cycling and non-cycling cells by comparison to a threshold γ s and to characterize the cell vulnerability towards apoptosis (programmed cell death). Phases 1, 2, 3 correspond respectively to ranges in the values of (a, γ) in the following open sets of the age-maturity plane: (see Figure 2):

Ω 1,k =]ka 2 , ka 2 + a 1 [×]0, γ s [, Ω 2,k =]ka 2 + a 1 , (k + 1)a 2 [×]0, γ s [, Ω 3,k =]0, (k + 1)a 2 [×]γ s , γ max [
for integers k = 0, 1, . . .. We will also use the notation

Q j,k = Ω j,k ×]0, T [ for j = 1, 2, 3.
The cell population in a follicle, f , is represented by cell density functions, φ k f,j (a, γ, t), defined on each cellular phase Q j,k , j = 1, 2, 3, k = 0, 1, . . . as solutions of the following conservation laws [START_REF] Echenim | Multi-scale modeling of the follicle selection process in the ovary[END_REF]: In phases 1 and 3, both a global control, U , and a local control, u f , act on the velocities of aging (g f function, locally controlled in phase 1) and maturation (h f function, locally controlled in phases 1 and 3) as well as on the loss term (λ apoptosis rate, globally controlled in phases 1 and 3). Phase 2 is uncontrolled and corresponds to completion of mitosis after a pure delay in age, a 2 -a 1 (no cell can leave the cycle here). We have the following expressions for those functions, for k = 0, 1, . . . (all parameters are real positive numbers defined in Table 1):

∂φ k f,j ∂t + ∂g f (u f )φ k f,j ∂a + ∂h f (γ, u f )φ k f,j ∂γ = -λ(γ, U )φ k f,j in Q j,k , j = 1, 2, 3, k = 0, 1, . . . (1)
g f (γ, u f ) = τ gf (1 -g 1 ω s-(γ)(1 -u f )) in Ω j,k , j = 1, 3 (2)
h f (γ, u f ) = τ hf (-γ 2 + (c 1 γ + c 2 )(1 -ω s (γ) exp(-u f /ū))) in Ω j,k , j = 1, 3 (3)
λ(γ, U ) = ω λ (γ)(1 -U ) in Ω j,k , j = 1, 3 (4)
g f (u f ) = τ gf , h f (γ, u f ) = λ(γ, U ) = 0 in Ω 2,k (5)
The ω s-and ω s+ functions are smooth switches between 0 and 1 in a strip defined for all a, 0 < γ s-< γ < γ s+ containing γ s (we use the notations |x| + = max(x, 0) and |x| -= max(-x, 0)):

ω s-(γ) = exp - |γ -γ s-| + |γ -γ s | , ω s+ (γ) = exp - |γ -γ s+ | - |γ -γ s | , ω s (γ) = ω s-(γ) + ω s+ (γ).
The factor ω s-in (2), switches off the control of g f in phase 3: g f = τ gf in Ω 3,k , and switches it on for (a, γ) strictly in phase 1, where

g f (u f ) = τ gf (1 -g 1 + g 1 u f )) for γ ≤ γ s-.
In the same manner, ω s in (3), switches off the control of h f for γ = γ s , so that, in a neighbourhood of γ s , h f > 0 for all t and u f : the boundary γ = γ s is inward for Ω 3,∞ and outward for each Ω 1,k . In particular, a cell very close to exit the cycle (in phase 1) is committed to go into phase 3. The exiting flux is controlled by u f inside phase 1 before the cell maturity reaches γ s : the sign of h f can be changed by the control only below γ s-.

In complement to this local control (by u f ) of the exit flux from the cell cycle, a global control is exerted (by U , U ≤ 1) on the cell vulnerability towards apoptosis, in a zone surrounding γ s , through the function

ω λ (γ) = K exp -γ-γs γ 2 in (4).
Each boundary segment of each (rectangular) Ω j,k , j = 1, 3, k = 0, 1, . . . is either inward or outward, independently on time t or controls u f , U , so that we can consider the following boundary conditions (the use of the trace values and the wellposedness of this model will be discussed later):

Boundary conditions for φ k f,1 , t ∈]0, T [ :

g f (γ, u f )φ k f,1 (ka 2 , γ, t) = 2τ gf φ k-1 f,1 (ka 2 , γ, t), for k ≥ 1, 0, for k = 0 for γ ∈]0, γ s [. ( 6
)
φ k f,1 (a, 0, t) = 0, for a ∈]ka 2 , ka 2 + a 1 [ (7)
Boundary conditions for φ N f,3 , N ≥ 1, t ∈]0, T [ :

φ N f,3 (0, γ, t) = 0, for γ ∈]γ s , γ max [ (8) φ N f,3 (a, γ s , t) = φ k f,1 (a, γ s , t), for a ∈ [ka 2 , ka 2 + a 1 [ 0, for a ∈ [ka 2 + a 1 , (k + 1)a 2 [ , for k = 0, . . . N -1. (9)
Initial conditions:

φ k f,j (a, γ, 0) = φ k f 0 (a, γ)| Ω j,k , j = 1, 3. (10)
In [START_REF] Besson | Existence and uniqueness of solutions for linear conservation laws with velocity field in L ∞[END_REF], the delay a 2 -a 1 and the doubling of the flux due to mitosis have been explicitly taken into account, so that it is not necessary anymore to consider phase 2 and φ k f,2 in the model. The cell density of interest in a follicle, φ N f , is simply:

φ N f = φ k f,1 on Q 1,k , k = 0, . . . N -1, φ N f = φ N f,3 on Q 3,N
In practice, we will fix the value of N to a large enough integer, and use the notations φ f and a max = N a 2 . φ f can be computed by solving the series of problems:

1. P 1,k : (1), ( 6), ( 7), [START_REF] Bryson | Applied optimal control[END_REF] on Q 1,k , for k = 0, . . . N -1. (11) 2. P 3,N : (1), ( 8), ( 9), [START_REF] Bryson | Applied optimal control[END_REF] on Q 3,N , when the P 1,k are solved .

(

) 12
It is worth noticing that those problems can be solved in sequence when the traces of the solutions of P 1,k on their outward boundaries (left and upper segments) are well defined. This is a useful property for the efficiency of computations and the backward reachability technique. The feedback exerted by the ovaries on the secretion of the hormonal control FSH defines a closed-loop system (cf. Figure 3). The global control U can be interpreted as FSH plasmatic level. The local control u f represents intra-follicular bioavailable FSH levels. It is given as a proportion of U . Define the maturity operator M as:

M (ϕ)(t) = γmax 0 amax 0 γϕ(a, γ, t)dadγ (13)
The global maturities M (φ f ) on the follicular scale, and M ( f φ f ) on the ovarian scale, will be used to define the two-scale feedback control:

U = S τ (M ( f φ f )) + U 0 u f = b f τ f (M (φ f )).U (14)
where

S(µ) = U s + 1 -U s 1 + exp (c(µ -m)) , S τ (µ)(t) = S(µ(t -τ )) b f (µ) = b 1 + 1 -b 1 1 + exp(b 2 (b 3 -µ)) , b f τ f (µ)(t) = b f (µ(t -τ f ))
The decreasing sigmoid function S accounts for the ovarian negative feedback on FSH ; U 0 (t) is a potential exogenous entry. The increasing sigmoid function b f makes u f increase with the maturity of follicle f , up to its reaching FSH plasmatic values. The delay τ is introduced to take into account the time needed for FSH signal to reach the ovaries; and the local delay τ f , that can be neglected in practice (τ f ≪ τ ) is mainly here to facilitate the resolution on finite time intervals. The delayed versions of S and b f are denoted S τ and b f τ f respectively. The parameters are chosen such that U and b f are normalized: 0 < b 1 , U s < 1, so that: Ovulation is triggered when estradiol levels reach a threshold value M s . As estradiol secretion is related to maturity (see [START_REF] Echenim | Multi-scale modeling of the follicle selection process in the ovary[END_REF]), the ovulation time T s is defined as:

0 ≤ u f ≤ U ≤ 1 + U 0 . (15
T s = min{T | M ( f φ f )(T ) = M s } (16)
The follicles are then sorted according to their individual maturity. The ovulatory follicles are those whose maturity at time T has overpassed a threshold M s1 such as M s1 ≤ M s . The ovulation rate is computed as:

N s,s1 = Card{f | M (φ f )(T s ) ≥ M s1 } (17)
For instance, Figure 4 shows, on the age-maturity domain, the cell density of either an ovulatory follicle (left panel) or an atretic one (right panel) at ovulation time. The cell cycle is implemented on a periodic domain, where the age a is reset after the cells go through mitosis. The granulosa cell repartition in the ovulatory follicle is characterized by a roughly older age range, and a higher maturity and cell density ranges than in the atretic follicle.

Wellposedness of the model

We have obtained a compact expression of the model by using a density representation of the cell population instead of a discrete model with a huge number of cells. Even if reasoning about controlled trajectories of individual cells will be useful in the sequel, the global numerical simulations are based on this compact form. It is then useful to discuss the wellposedness of the model in order to check, to some extent, its consistency. Due to the delays in the feedback control laws, we can consider the u f and U control terms as given functions of time on each interval of size τ f . The velocities and loss term are then given functions g f (γ, t), h f (γ, t) and λ(γ, t). In this case, we show below that there exists a unique solution to the series of problems [START_REF] Mitchell | A Toolbox of Hamilton-Jacobi Solvers for Analysis of Nondeterministic Continuous and Hybrid Systems[END_REF], [START_REF] Mitchell | A toolbox of level set methods[END_REF]. In each domain Q 1,k , Q 3,N , the velocity field (g f , h f ) and the loss rate λ can be considered, without loss of expressiveness for the model, as smooth functions. Yet, initial and boundary values may be irregular, so that we need to check the existence of weak solutions. As mentioned below, it is sufficient to solve the sequence of P 1,k and then P 3,N to define φ f . Since those problems have a common structure, we just need to study the solution of P 1,k and its traces on the outward boundary used to define both P 1,k+1 and P 3,N .

We first use a change of variables to transform (1) into a conservation equation without a loss term. Let

φ k f,1 = φ k f,1 exp (l(γ, t)) defined on Q 1,k . Let ∂Q + 1
,k and ∂Q - 1,k denote the inward and outward boundaries of the domain Q 1,k , and let l(γ, t) such as

∂l(γ, t) ∂t + h f (γ, t) ∂l(γ, t) ∂γ = -λ(γ, t), l(γ, 0) = 0 (18)
Here, we can assume that h f and λ (extended to all γ ∈ R) are C 1 functions, so that, according to standard theory, there exists a unique generalized solution l of the linear transport equation (18), that is Lipschitz for all t > 0. Now φ

k f,1 verifies ∂ φ k f,1 ∂t + div(v f φ k f,1 ) = 0 (19)
where v f is the (smooth) velocity vector:

v f = g f (γ, t) h f (γ, t) .
Non homogeneous boundary conditions on ∂Q + 1,k , similar to (6), ( 7) complete (19). To solve such a problem, with L ∞ or L 1 boundary and initial values, extensions to the classical Kruzkov entropy solutions ( [START_REF] Kruzkov | First order quasilinear equations in several independent variables[END_REF]) to boundary value problems have been proposed (see e.g. [START_REF] Szepessy | Scalar conservation laws with boundary conditions and rough data measure solutions[END_REF] and [START_REF] Porretta | L 1 solutions to first order hyperbolic equations in bounded domains[END_REF] where renormalized entropy solutions are defined for this type of problem). An additional problem faced here consists in defining the trace of the solution on the outward boundary ∂Q - 1,k , in order to solve the next problems P 1,k+1 . It is not clear whether this trace can be defined here for such weak solutions. As

v f ∈ (L ∞ (Q 1,k )) 2 and div(v f ) ∈ L ∞ (Q 1,k
), an alternative is to use the results of [START_REF] Besson | Existence and uniqueness of solutions for linear conservation laws with velocity field in L ∞[END_REF], to define a unique solution φk f ∈ L 2 (Q 1,k ) with well defined trace in L 2 (∂Q - 1,k ) as soon as the boundary data are in L 2 (∂Q + 1,k ). A drawback of this alternative is that, in our case, it is not known whether those solutions coincide with classical weak solutions.

Characteristic equations

The model deals with controlled partial differential equations. Up to now, no convincing control strategy for such velocity controlled conservation laws with integro-differential control terms is available from the literature. To tackle the control problems, we focus on the actions of the control terms on the characteristics of the conservation laws. As the solution of the conservation laws is L 2 , we can derive the characteristic equations. We now deal with the following ordinary differential equations [START_REF] Evans | Partial differential equations[END_REF]:

     ȧc = g f (γ c , u f ) γc = h f (γ c , u f ) φ f c = -λ(γ c , U ) + ∂h f (γc,u f ) ∂γc φ f c c = 1 . . . n (20)
where g f , h f and λ are defined in (2),( 3),( 4), [START_REF] Porretta | L 1 solutions to first order hyperbolic equations in bounded domains[END_REF] and n is the number of characteristics.

The local control u f acts on the position of the characteristics on the (a c , γ c ) spatial domain at a given time, and U on the loss term λ.

The transfer conditions between the cellular phases are continuous on a c and γ c and discontinuous on φ f c .

From phase 1 to phase 2 (crossing the right boundary of some Ω 1,k ), the flux continuity condition reads:

-for T k such as:

a c (T k ) = a 1 + ka 2 , 0 ≤ γ c (T i ) ≤ γ s , for some k ∈ N, φ f c (T + k ) = 1 -g 1 ω s-(γ)(1 -u f (T - k )) φ f c (T - k ) (21)
From phase 2 to phase 1 (crossing the left boundary of some Ω 1,k ), after mitosis, the flux doubling condition reads:

-for T k such as:

a c (T k ) = ka 2 , 0 ≤ γ c (T i ) ≤ γ s for some k ∈ N, 1 -g 1 ω s-(γ)(1 -u f (T + k )) φ f c (T + k ) = 2φ f c (T - k ) (22)
Ovulatory follicles are distinguished from the atretic ones by their characteristic curves reaching a well defined zone of the (a c , γ c , φ f c ) space at a given time, as we can see on Figure 4. In the next section, we study to which extent characteristics curves can reach the target sets when subject to appropriate open-loop control laws.

3 Backwards reachable sets 3.1 Link with Hamilton-Jacobi-Bellman equations

We assume that a follicle ovulates (respectively undergoes atresia) at time

T if: ∀c = 1 . . . n, (a c (T ), γ c (T ), φ f c (T )) ∈ M o ( resp. ∈ M a )
where M o and M a are the targets for respectively ovulation and atresia.

We focus here on the backwards problem, and we first introduce the elements of its solution, the solvability tube [START_REF] Kurzhanski | Dynamic optimization for reachability problems[END_REF] (or the backward reachable set [START_REF] Mitchell | A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games[END_REF]) associated with M o (resp. M a ), i.e. the states from which it is possible to reach M o (resp. M a ), given admissible controls U and u f . From a physiological viewpoint, it corresponds to the initial conditions compatible with ovulation (resp. atresia) and subject to correct control laws amongst the set of admissible controls.

We recall briefly how such backwards reachability problems can be solved in the framework of optimal control theory. Let us consider the following nonlinear continuous controlled dynamics, where f is a Lipschitz continuous function, so that there is a unique continuous solution for any measurable u(t):

   ẋ = f (t, x, u), t ≥ τ x(τ ) ∈ R n u ∈ U ⊂ R m (23)
Given a closed target set M ∈ R n and two times τ and t 1 the solvability set W (τ, t 1 , M) is the set of states x ∈ R n such that there exist control functions u with values in the compact set U that steer system (23) from the state x(τ ) = x to x(t 1 ) ∈ M [START_REF] Kurzhanski | Dynamic optimization for reachability problems[END_REF]. W (τ, t 1 , M) can be computed by solving the following optimization problem where d 2 (x, X) = min{||x -z|| 2 | z ∈ X} is the square of the distance d(x, X) from point x to set X:

min u∈U {d 2 (x(t 1 ), M)|x(τ ) = x} (24)
The solution of this optimal control problem can be found by solving the following backward Hamilton-Jacobi-Bellman (HJB) equations, where

H f (D x V, x) = min u∈U {(D x V ) T .f (t, x, u)} is the Hamiltonian [10]: V t + H f (D x V, x) = 0, t ≤ t 1 , V (t 1 , x) = d 2 (x(t 1 ), M) (25)
Then the following property is true [START_REF] Kurzhanski | Dynamic optimization for reachability problems[END_REF]:

W (τ, t 1 , M) = {x|V (τ, x) ≤ 0} (26)
The solvability tube is then the set-valued function of the starting time t for given t 1 and M:

W[t 1 , M] : t → W[t 1 , M](t) = W (t, t 1 , M), τ ≤ t ≤ t 1 . (27)
Any initial time t is thus associated with a delay of exactly (t 1 -t) to reach the target. In contrast, the backward reachable set is the set of states from which it is possible to reach the target in at most (t 1 -t) [START_REF] Mitchell | A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games[END_REF]:

G(t ; t 1 , M) = t≤s≤t1 W (t, s, M), t ≤ t 1 . (28)
This set can also be characterized in the same manner as in expression (26). It is the zero sublevel set of the V function solution of the HJB equation:

V t + min(0, H f (D x V, x)) = 0, t ≤ t 1 , V (t 1 , x) = d 2 (x(t 1 ), M) (29) G(t ; t 1 , M) = {x ∈ R n |V (t, x) ≤ 0} (30)
In the next section, we use the numerical methods based on equations ( 29), (30) to compute the G set [START_REF] Mitchell | A Toolbox of Hamilton-Jacobi Solvers for Analysis of Nondeterministic Continuous and Hybrid Systems[END_REF][START_REF] Mitchell | A toolbox of level set methods[END_REF]. This choice is convenient since it seems easier to compute the G set than the W set-valued function, but it amounts to deal with desynchronized final times. Indeed, in the W[t 1 , M] tube, every states considered at the same starting time reach the M target at the same final time t 1 , whereas in G(t ; t 1 , M) states reaching the M target at different final times (provided that they are not later than t 1 ) are fused. To ensure that every states reaching the target remain inside until t 1 , the control function has to be extended.

Application to granulosa cells

The follicle as a controlled dynamical system

We represent each follicle by n cells following characteristic curves of (1) defined by three state variables: age, maturity and cell density as in (20). The model of the follicle is then of the form (23) with 3n state variables:

x(t) = (a 1 (t), ..., a n (t); γ 1 (t), ..., γ n (t); φ f 1 (t), ..., φ f n (t)) T and for each cellular phase i = 1, 2, 3, the dynamics is defined by the following f i :

f i =                  g f (γ 1 , u f ) . . . g f (γ n , u f ) h f (γ 1 , u f ) . . . h f (γ n , u f ) K(γ 1 , u f , U )φ f 1 . . . K(γ n , u f , U )φ f n                  (31)
where g f and h f are defined in (2), (3), (5) and K(γ c , u f , U ) is defined by Phases 1 and 3:

∀k, N ∈ N, ∀(a c , γ c ) ∈ Ω 1,k ∪ Ω 3,N K(γ c , u f , U ) = -λ(γ c , U ) + ∂h f (γ c , u f ) ∂γ c (32) Phase 2: ∀k ∈ N, ∀(a c , γ c ) ∈ Ω 2,k K(γ c , u f , U ) = ln(2)δ(a c )φ f c , with (k+1)a2 ka2+a1 δ(a)da = 1 (33)
The expression of K(γ c , u f , U ) in phase 2 results from the regularization by a smooth positive function δ(a) of the mitosis at the transition (22) between phases 2 and 1, that will be used in the numerical treatment of the equations where the transmission boundaries are included in the computational domain. The precise choice for δ(a) does not matter since its only role is to guarantee that the density of cells entering phase 2 at time ka 1 has doubled at time ka 1 + a 2 . Since a 2 -a 1 = 1, the simplest choice is δ(a) = 1 so that we have

δ(a c ) = 1, resulting in K(γ c , u f , U ) = ln(2)φ f c .
According to physiological knowledge, we can consider that there exists a maximal age, a max , beyond which all cells become senescent and die. This assumption introduces the notion of a maximal life time T max such that a(T max ) = a max (T max exists because ȧ > g 2 > 0). Hence, there is also a maximal maturity γ max and a maximal cell density φ f max attainable for t ∈ [0, T max ]. The f i are Lipschitz continuous, so that there is a unique solution x within each phase.

We now state some properties of this controlled system and the corresponding reachable targets in the age-maturity plane. We first describe the control action on the orientation of the velocity field:

In the interior of phases 1 & 3:

∀u f ∈ U, g h ≥ g 1 τ gf ≡ 0.5τ gf (34)    h f (γ, u f ) ≥ 0, for u f ≥ u * f and 0 ≤ γ ≤ γ + (1/ū) h f (γ, u f ) ≤ 0, for u f ≤ u * f and 0 ≤ γ ≤ γ + (1/ū) h f (γ, u f ) ≤ 0, for all u f and γ + (U/ū) ≤ γ ≤ γ max (35)
In the interior of phase 2:

∀u f ∈ U, g h = 1, h f = 0. ( 36
)
with:

u * f (γ) =    ū ln c 1 γ + c 2 c 1 γ + c 2 -γ 2 for 0 ≤ γ ≤ γ + (1/ū) 1 for γ + (1/ū) ≤ γ ≤ γ max (37) γ ± (ν) = c 1 (ν) ± c 1 (ν) 2 + 4c 2 (ν) 2 , with c i (ν) = c i (1 -exp(-ν)), i = 1, 2. ( 38
)
We just sketch the proof (easy but tedious) of these properties. We can first remark that, in the interior of phases 1 and 3 excluding the strip defined by γ s± , the ω s± functions can only take as values 0 or 1. In this interior subset h f (γ, u f ) is defined as

h f (γ, u f ) = τ hf (-γ 2 + (c 1 γ + c 2 )(1 -ω s (γ) exp(-u f /ū))) = τ hf (γ + (u f /ū) -γ)(γ -γ -(u f /ū))
It is easy to check that h f has the sign of (γ + (u f /ū) -γ). Besides, γ + is an increasing function of u f , since

γ ′ + (ν) = c 1 γ + (ν) + c 2 c 1 (ν) 2 + 4c 2 (ν) exp(-ν) Hence 0 ≤ γ + (u f /ū) ≤ γ + (u max /ū) ≡ γ + (1/ū), so that we can solve h f (γ * , u f ) = 0 and define u * f (γ) such that γ + (u * f (γ)/ū) = γ for 0 ≤ γ ≤ γ + (1/ū). Finally, as exp - 1 ū is very small, γ + (1/ū) ≈ (1 -exp(-1/ū)) c 1 + c 2 1 + 4c 2 2 < γ max
This yields property (35). The remaining properties (34) and (36) are obvious.

Such results imply the following properties of the reachable targets:

1. As g h > 0, the ages of cells in the reachable targets set have to be larger than those of the cells in their initial conditions. Hence, the reachable targets are on the right of the cell initial positions in the age-maturity plane.

2. For a given global control U and maturities 0 ≤ γ ≤ γ + (U/ū), property (35) shows that it is possible, using a local feedback u f around the feedback u * f , to either increase, decrease or maintain the maturity.

In particular, if a target maturity is reached and U remains large enough, the cells remains inside the target. For such targets M, G(t ; t 1 , M) = W (t ; t 1 , M). It is also possible to reach a target and then leave it, e.g. if U is lowered. This set-equality depends then strongly on the control strategies and positions of the targets with respect to γ + (1/ū).

3. As U ≤ 1, maturities γ + (1/ū) ≤ γ ≤ γ max are not reachable.

Subsequently, we assume that:

the target ages are large enough (e.g. > a 1 + a 2 ) (39) the target maturities are below γ + (1/ū) (40) so that it is not impossible to reach them (but success is not guaranteed, the density conditions needing also to be met).

The Hamiltonian in each cellular phase

The HamiltonJacobi-Bellman (HJB) equations associated with the backwards reachable sets are solved once a target has been defined and the Hamiltonians in each phase have been minimized according to the control terms. We do not use here the feedforward term: U 0 = 0. Hence both the global and local FSH signals operate as best-case feedback controls subject to the explicit constraints (rewriting (15

)) 0 ≤ u f ≤ U ≤ 1.
We note u = (u f , U ) T the control and U the set of admissible controls. Such constraints take into account the following physiological knowledge:

(i) the level of FSH in the antral 1 fluid is at much as high as FSH plasmatic levels, which implies u f ≤ U ;

(ii) there is a basal secretion (i.e. not subject to ovarian feedback) of FSH by the pituitary gland, which implies U ≥ 0;

(iii) the balance between the secretion rate (or exogenous administration rate) and the clearance rate ensure the saturation of FSH plasmatic levels, which implies U ≤ 1. Two coupled biological questions underlie this reachability problem: (1) which patterns of exogenous FSH administration can be applied to target either ovulation or atresia and (2) how follicular vascularization should develop to remain compatible with those patterns. The answer to such questions might help improving the current FSH treatments, which suffer from several drawbacks such as ovarian overstimulation syndrome.

To reach a given target M from its current location in phase Ω i,k 2 , a cell (or the lineage of its daughter cells) has to reach some intermediate or final target in the closure of the Ω i,k domain. For instance, in the case of the ovulatory target

M = M o (M o ⊂ Ω 3,N ) W[t 1 , M o ]
following the solvability tube backward from Phase 3 leads to the boundary between Ω 3,k and Ω 1,k (k ≤ N ); the intersection of this boundary with the solvability tube defines a new target set M o,1,k , in Ω 1,k . Such intersections depend smoothly upon the width of the thin strip 0 < γ s-< γ < γ s+ containing these boundaries where the control has little effects. In the following, we suppose that γ s-= γ = γ s+ so that ω s- and ω s+ are replaced by perfect switch functions.

In each cellular phase i, in order to compute the backward reachable set G(t 1 -t ; t 1 , M i ), i.e. the set of states from which it is possible to reach the target M i in less than a given time t, we have to solve an HJB equation of the form (29) and use the definition (30):

V t + min(0, H f i (D x V, x)) = 0, t ≤ t 1 , V (t 1 , x) = d 2 (x(t 1 ), M i ) (41) G(t ; t 1 , M) = {x ∈ R n |V (t, x) ≤ 0} (42)
where

H f i (p, x) = min u∈U (p T .f i (x, u)), with p = D x V = (p a1 , ..., p an , p γ1 , ..., p γn , p φ f 1 , ..., p φ f n ) T .
1 ovarian follicles are spheroidal structures hollowed by a cavity called the antrum 2 where i stand for the phase index and k for the number of cell generation elapsed since initial time

Phase 2

This phase is not controlled so that the expression of the Hamiltonian is straightforward using (33):

H f 2 = τ gf n c=1 p ac + ln(2) n c=1 p φ f c δ(a c )φ f c
This analytical minimization enables to write the expression of the Hamiltonian in each cellular phase.

Numerical results

Implementation: level set methods

Two main approaches are used to find reachable sets: the Eulerian approach approximates the solution values on a fixed grid, whereas the Lagrangian approach tracks the trajectories of the dynamics. In case of backwards reachable sets, Eulerian methods are the most appropriate, as they enable to correctly compute the solution beyond shocks [START_REF] Mitchell | Application of level set methods to control and reachability problems in continuous and hybrid systems[END_REF].

Amongst the various Eulerian methods available to compute the exact reachable sets, we used "level set methods" to solve Eq.( 25). They enable to simulate the motion of dynamic surfaces, such as the zero-level set of HJB equations [START_REF] Mitchell | A toolbox of level set methods[END_REF], with a very high accuracy (about a tenth of the spacing between the grid points [START_REF] Mitchell | Application of level set methods to control and reachability problems in continuous and hybrid systems[END_REF]). Such level set algorithms are implemented in the "Toolbox of Level Set Methods"3 , which can be used to solve, among others, equations of the form:

D t V (x, t) + H(x, D x V ) = 0 with the constraints D t V (x, t) ≥ 0 or V (x, t) ≥ Ṽ (x, t)
The time derivative, D t V , is approximated by a Runge-Kutta scheme, with customizable accuracy. The spatial derivative, D x V , is approximated with an upwind finite difference scheme. The Hamiltonian, H(x, p), is approximated with a Lax-Friedrichs scheme:

Ĥ(x, p + , p -) ≡ H(x, p + + p - 2 ) - 1 2 α T (p + -p -)
where p + and p -are respectively the right and left approximations of p, and α is an artificial dissipation coefficient added to avoid oscillations in the solution.

Although the toolbox is designed for solving initial value problems, it is possible, in case of autonomous systems, to solve final value problems by reversing the time in system (23) [START_REF] Mitchell | A toolbox of level set methods[END_REF].

Simulation results

The initial problem is a 3 × n dimensional problem, where n is the number of characteristics considered. This problem is not tractable from a numerical ground. Indeed, for n ≥ 2 and a correct number of grid points in each dimension (≃ 100) we cope with too many points (100 3n ), and face the dimensionality limits of the algorithms. We thus computed the backwards reachable sets for one characteristic, (a 1 , γ 1 , φ f 1 ), corresponding to an ovulatory or atretic trajectory. It is a simplified case where all granulosa cells in a follicle are assumed to be synchronized in age and maturity.

In the case of ovulatory follicles, the boundaries of the target set M o are chosen to coincide with the ranges of age, maturity, and cell density reached by the characteristics of an ovulatory follicle at ovulation time, (see Figure 4, left).

The target set M o for ovulation is thus defined as: This is in agreement with the cell density dynamics in phase 3.

M o = {(a 1 , γ 1 , φ f 1 ) ∈ [
The backwards reachable set of M o nearly contains the whole spatial domain. It is an overapproximation of the initial states that allow to reach the target set. Indeed, it does not only contain the trajectories of the characteristics of the ovulatory follicle, but also those starting from initial conditions that do not arise on a physiological ground. Basing on biological knowledge, we can thus select from the overapproximated set, the physiologically-meaningful subset: we assume that all cells are initially within the cell cycle, at their first generation, so that the subset of admissible initial conditions is defined by {(a

1 (t 0 ), γ 1 (t 0 )) ∈ [0, a 2 ] × [0, γ s ]}.
This subset is below the dashed segment on the left panel of Figure 5. The target set M a for atresia is defined by the ranges reached by the characteristics of atretic follicles:

M a = {(a 1 , γ 1 , φ f 1 ) ∈ [8, 10] × [3, 4] × [2, 4]} (48)
Figure 6 illustrates all the trajectories that fall into the target set for an atretic follicle. The left panel illustrates the projection, on the age-maturity plane, of the changes in the backwards reachable set. Some states in the target set are initially in the cell cycle, leading to the "staircase" pattern described below. The right panel illustrates the projection, on the age-cell density plane, of the temporal changes in the backwards reachable set. The set area increases, and at time t = 11, it roughly contains the whole domain. Once again, the set of initial states leading to atresia is overapproximated, and we can apply the same constraints as in the ovulatory case on the initial conditions, that are again below the dashed segment on Figure 6. The initial conditions leading to ovulation or atresia, subject to physiological constraints, are roughly superimposed, and contain the whole subset of admissible initial conditions. Thus we cannot distinguish a priori ovulatory trajectories from atretic ones. This result agrees with physiological knowledge, as there is no predestination in follicular fate [START_REF] Gougeon | Regulation of ovarian follicular development in primates: facts and hypotheses[END_REF].

Although the initial conditions do not enable to determine whether a follicle is ovulatory or atretic, there is a moment when trajectories leading to either ovulation or atresia diverge. Such a moment corresponds to the physiological time of selection, when the trajectories of the few ovulatory follicles separate from the atretic ones.

Discussion and perspectives

We have characterized the backwards reachable sets for both ovulatory and atretic follicles. This control problem has been solved by studying the characteristics of the conservation laws describing the dynamics of the structured granulosa cell populations in ovarian follicles. The HJB equations representing the backwards reachable sets have been solved numerically for one characteristic curve, and allow to characterize the set of initial conditions that lead a characteristic into the target set either for ovulation or atresia.

We have used the backwards reachable sets theory for a problem that is not entirely continuous, as there are discontinuities at the transitions between each cellular phase. Hence, we studied each phase separately, assuming that the transitions at the boundaries were continuous. Actually, we deal with a hybrid problem, as each cellular phase has its proper dynamics. As far as time dependent HJB equations are concerned, elements of a theory for hybrid systems are established, especially concerning "reach-avoid" sets [START_REF] Tomlin | Computational techniques for the verification of hybrid systems[END_REF], but it is not complete yet, and not directly applicable to our problem. Viability theory also considers reachable sets for hybrid systems [START_REF] Aubin | Impulse differential inclusions and hybrid systems: a viability approach[END_REF], but the algorithms used are not as accurate as level set methods, since they use a different representation of reachable sets [START_REF] Mitchell | Application of level set methods to control and reachability problems in continuous and hybrid systems[END_REF]. Yet the numerical results we have obtained with the level set methods seem correct since the continuous part of the system is solved with highly accurate algorithms and the hybrid part is handled via continuity conditions and substantiated by the trace results obtained on the conservation laws.

We have solved a simplified numerical problem dealing with a single characteristic of a follicle, as the simulation tool needs too huge memory to solve higher dimensional problems. The backwards reachable sets obtained show that it is possible to find a correct control law to steer a large set of initial conditions into the target set, and they give information upon the maximum duration needed to reach the target.

If the simulation of more than one characteristic were possible, we expect that the backwards reachable sets would be smaller. We even speculate that, from a physiological viewpoint, a larger set in case of a single growing follicle ensures the ovulation success, while a smaller set in case of several simultaneously growing follicles limits the ovulation rate and underlies the selection process in part. Yet this assumption needs improvement of the numerical implementation to be tested.

We have studied by now the open-loop problem, for the ovulation of one follicle. This situation is close to therapeutic situations of ovarian stimulations, so that the designed control laws may be useful in improving therapeutic schemes. Yet, to understand ovulation control in depth, we will have to tackle multi-scale, closed-loop reachability problems of selection, which is the matter of future work.

FSHFigure 1 :

 1 Figure 1: Cell flow chart. The cell cycle consists of the cyclic G1-SM-G1 pathway. When it enters G1 from SM, a mother cell gives birth to two daughter cells. Cell differentiation corresponds to the one-way flow from G1 into D. Entry into apoptosis arises from the G1 and D phases. The G1 and D phases are under the control of FSH signal.

Figure 2 :

 2 Figure 2: Cellular phases on the age-maturity plane. The cell cycle is represented by the pathway Ω 1,k , Ω 2,k ... and the differentiation phase by Ω 3,N .

Figure 3 :

 3 Figure 3: Ovarian system: the multi-scale model in closed-loop. The global feedback control U is modulated on the follicle scale (F ol i), by the follicular maturity, M (φ i ), to define the local control u f i .

Figure 4 :

 4 Figure 4: Cell density of an ovulatory follicle (left) and an atretic follicle (right). The horizontal axis represents the cell age, the vertical axis the cell maturity, and the colorbar indicates the cell density value (number of cells per elementary age × maturity volume). a 2 is the cell cycle duration (2 age units).

Figure 5 :

 5 Figure 5: Solvability tubes for ovulatory follicles. The left panel shows the projection of the 3D set on the (a 1 -γ 1 ) plane, and the right one on the (a 1 -φ f 1 ) plane.

Figure 6 :

 6 Figure 6: Solvability tubes for atretic follicles. The left panel shows the projection of the 3D set on the (a 1 -γ 1 ) plane, and the right one on the (a 1 -φ f 1 ) plane.

Table 1 :

 1 ) Main model functions and parameters

	Functions and parameters Definition

  Figure5illustrates the trajectories compatible with ovulation. The left panel illustrates the projection, on the age-maturity plane, of the changes in the backwards reachable set. Time t = 0 (first subplot) is ovulation time, and we can see the rectangle {[START_REF] Bryson | Applied optimal control[END_REF][START_REF] Mitchell | A toolbox of level set methods[END_REF]×[START_REF] Bryson | Applied optimal control[END_REF][START_REF] Mitchell | A Toolbox of Hamilton-Jacobi Solvers for Analysis of Nondeterministic Continuous and Hybrid Systems[END_REF]} representing the projection of the target set M o (the age range is slightly different from that of Figure4due to difference in the dimensioning of the unrolled domain compared to the periodic domain). At time t = 4, the set area has increased, including all states that can join the target set in less than 4 time units (2 time units roughly correspond to 1 cell cycle duration). The backwards reachable set at time t = 11 represents the set of states that can reach the target set for ovulation in at most 11 time units. We chose to stop the simulation at this time as it contains the initial conditions used for the simulation represented in Figure4. The "staircase" shape of the backwards reachable set is due to the dynamics in the cell cycle. The maturity increases with age in phase 1, while its remains unchanged in phase 2, yielding to the flat part of the staircase. An instance of the flat part can be seen at time t = 11 between a 1 = 7 and a 1 = 8. The right panel illustrates the projection, on the age-cell density plane, of the changes in the backwards reachable set. We can again see the projection of the target set M o at time t = 0, as the rectangle {[START_REF] Bryson | Applied optimal control[END_REF][START_REF] Mitchell | A Toolbox of Hamilton-Jacobi Solvers for Analysis of Nondeterministic Continuous and Hybrid Systems[END_REF] ×[START_REF] Szepessy | Scalar conservation laws with boundary conditions and rough data measure solutions[END_REF][START_REF] Besson | Existence and uniqueness of solutions for linear conservation laws with velocity field in L ∞[END_REF]}. At time t = 11, a large part of the age-cell density domain is compatible with reaching the target set. From age a 1 = 9 onward, the cell density variable increases with age to reach the target set.

	10, 12] × [10, 11] × [4, 6]}	(47)

† Unite de Recherche INRIA Rocquencourt, Domaine de Voluceau, Rocquencourt BP 105, 78153 Le Chesnay Cedex, France tel: +33 1 39 63 57 41 fax: +33 1 39 63 57 86

http://www.cs.ubc.ca/˜mitchell/ToolboxLS/

If the targets fulfil the conditions (39) and (40) discussed above, we can conclude that: i) Starting from cell ages "younger" than the target ages, for all u ∈ U, the cells get closer in age to target: For all i, p ai .g f (γ i , u f ) ≤ 0, as p ai ≤ 0 and g f > 0. ii) For the particular choice U = 1 and u f = u * f , we have h f (γ i , u f ) = 0, so that: For all i, p γi .h f (γ i , u * f ) = 0. In particular, p T .f i (x, u * f )) ≤ 0, so that H f i (p, x) ≤ 0. We have thus shown that in our case, solving equation (41) amounts to solve the following simpler HJB equation to compute G:

We now compute the Hamiltonian H f i (p, x) in each phase.

Phases 1 & 3

We first compute p T .f i (x, u), for i = 1, 3.

which can be rewritten as follows:

with (using equations ( 2), ( 3) and ( 32))

where

Now we minimize this expression with respect to u ∈ U.

It is equivalent to minimize first over U with u f ≤ U ≤ 1 and then over u f with 0 ≤ u f ≤ 1 (we drop the arguments x, p for simplicity):

When C < 0, we have

We have then to solve min