Alexander Lorz
email: alexander.lorz@upmc.fr

Tommaso Lorenzi
email: tommaso.lorenzi@upmc.fr

Jean Clairambault
email: jean.clairambault@inria.fr

Benoît Alexandre Escargueil
email: alexandre.escargueil@inserm.fr

Benoît Perthame
email: benoit.perthame@upmc.fr

Modeling the effects of

Introduction

The levels of abiotic factors are not uniform throughout solid tumors [START_REF] Olive | Inhibition of Hedgehog signaling enhances delivery of chemotherapy in a mouse model of pancreatic cancer[END_REF][START_REF] Rainey | Adaptive radiation in a heterogeneous environment[END_REF][START_REF] Vaupel | Blood supply, oxygenation status and metabolic micromilieu of breast cancers: characterization and therapeutic relevance[END_REF][START_REF] Vaupel | Blood flow, oxygen and nutrient supply, and metabolic microenvironment of human tumors: a review[END_REF]. This leads to the creation of distinct niches differentiated by the local environment, which provide ecological opportunities for diversification (e.g., cells characterized by different proliferation or drug-resistant phenotypes can be selected depending on the local concentrations of nutrients and anti-cancer drugs), and paves the way for the emergence of intra-tumor heterogeneity [START_REF] De Bruin | Intra-tumor heterogeneity: lessons from microbial evolution and clinical implications[END_REF][START_REF] Gerlinger | Intratumor heterogeneity and branched evolution revealed by multiregion sequencing[END_REF]. In this line, histopathological evidence supports the idea that the emergence of resistance to anti-cancer therapies can be considered as a process of Darwinian-like selection in tumor cell populations [START_REF] Gerlinger | Intratumor heterogeneity and branched evolution revealed by multiregion sequencing[END_REF][START_REF] Merlo | Cancer as an evolutionary and ecological process[END_REF]. In fact, malignant clones with heterogeneous genetic/epigenetic expressions leading to different phenotypes (e.g., epithelial vs mesenchymal, with the same genetic material [START_REF] Weinberg | The Biology of Cancer[END_REF]) can be seen as competing for space and resources (i.e., oxygen, glucose or other nutrients) within the environment defined by the surrounding tissues, together with the selective pressure exerted by therapeutic actions.

Cytotoxic drugs, such as DNA damaging agents and antimetabolites, are the most widely used chemical weapons in the fight against cancer. However, these drugs usually cause unwanted toxic side effects in the patients' organisms, since they are designed to kill cancer cells and are seldom specific toward tumor cells. Furthermore, they tend to kill strongly proliferative clones, usually considered as made of the most drug-sensitive cells [START_REF] Mitchison | The proliferation rate paradox in antimitotic chemotherapy[END_REF], thus selecting for cells resistant to therapies [START_REF] Goldie | Drug resistance in cancer: mechanisms and models[END_REF][START_REF] Gottesman | Mechanisms of cancer drug resistance[END_REF][START_REF] Komarova | Drug resistance in cancer: Principles of emergence and prevention[END_REF][START_REF] Szakács | Targeting multidrug resistance in cancer[END_REF]. These are the two major obstacles -toxic side effects and emergence of resistant clones -encountered in the clinic when using cytotoxic agents in treating tumors. This situation calls for therapy optimization, that is, identification of drug doses and design of optimal delivery schedules in multi-drug combinations, allowing for an effective control of cancer growth. As regards multi-drug combinations, a trend in modern treatment of cancers leads to combining cytotoxic and cytostatic drugs. Cytostatic drugs slow down proliferation of cancer cells by blocking growth factor receptors, or downstream intracellular pathways involved in proliferation (e.g., tyrosine kinase inhibitors) and have lower toxicity for healthy cells. This mode of therapy combination allows to attain the twofold goal of reducing toxicity and holding in check the multiplication of resistant clones [START_REF] Silva | A theoretical quantitative model for evolution of cancer chemotherapy resistance[END_REF][START_REF] Tomasetti | Drug resistance always depends on the turnover rate[END_REF][START_REF] Tomasetti | An elementary approach to modeling drug resistance in cancer[END_REF][START_REF] Tsai | The molecular basis of targeting protein kinases in cancer therapeutics[END_REF][START_REF] Weinstein | Oncogene addiction[END_REF].

As far as drug delivery schedules are concerned, it has been suggested that, among others, infusion protocols based on bang-bang control (i.e., those protocols in which drug delivery is alternatively switched on and off over time) can allow an effective control of tumor size [START_REF] Łędżewicz | Optimal bang-bang controls for a two-compartment model in cancer chemotherapy[END_REF]. We will consider here the case where tumor cells are exposed to rectangle-wave infusions of cytotoxic and/or cytostatic drugs at constant concentrations.

Motivated by these considerations, and focusing on an in vitro cancer cell population as reference system (i.e., not taking into account unwanted toxicity to healthy cells, which is a theme we had explored in a previous paper [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF]), we propose a space and phenotype structured population model for selection dynamics under the effects of cytotoxic and/or cytostatic drugs. The model we design includes birth and death processes of cancer cells. Furthermore, it also takes explicitly into account the dynamics of resources and anti-cancer drugs as well as their interactions with the cell population under treatment.

Cancer cells are assumed to be organized in a non-vascularized and radially symmetric spheroid. The assumption of radial symmetry is consistent with experimental evidence on tumor spheroids of small size [START_REF] Yu | Holographic optical coherence imaging of rat osteogenic sarcoma tumor spheroids[END_REF]. As a simplification, we state our model in one space dimension [START_REF] Byrne | Continuum models of avascular tumor growth[END_REF][START_REF] Jackson | A mechanical model of tumor encapsulation and transcapsular spread[END_REF]. In spite of this apparently simplified assumption, we do believe our model relevant to study small small avascular tumors and micrometastases [START_REF] Carver | Multicellular tumor spheroids as a model for assessing delivery of oligonucleotides in three dimensions[END_REF][START_REF] Sutherland | Cell and environment interactions in tumor microregions: the multicell spheroid model[END_REF], which are believed to be the target of adjuvant chemotherapies, but also to tackle, as an initial step, highly complex macrotumor models targeted by metastatic chemotherapy [START_REF] Beyer | Coadministration of epithelial junction opener JO-1 improves the efficacy and safety of chemotherapeutic drugs[END_REF][START_REF] Beyer | Epithelial junction opener JO-1 improves monoclonal antibody therapy of cancer[END_REF].

Cells are structured as a population by two non-negative real variables x ∈ [0; 1] and r ∈ [0; 1] standing, respectively, for the normalized expression level of a cytotoxic resistant phenotype and for the normalized linear distance from the center of the spheroid. This implies that we do not consider a cell to be necessarily either totally sensitive or totally resistant to a given drug; we rather introduce a continuous structuring variable describing resistance between 0 (highly sensitive) and 1 (highly resistant).

It should be noted that, compared to the model proposed in [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF], the present one is able to mimic the simultaneous selection of several traits (i.e., the rise of phenotypic polymorphism) within the cancer cell population, which provides the basis for intra-tumor heterogeneity. The additional spatial structure variable r, together with the diffusion along the r-axis of nutrients and therapeutic drugs, are the key ingredients of this model that enable the emergence of such a heterogeneous scenario, which is close to the ones observed in biological experimentations [START_REF] Busch | Fluence rate-dependent intratumor heterogeneity in physiologic and cytotoxic responses to Photofrin photodynamic therapy[END_REF][START_REF] Gerlinger | Intratumor heterogeneity and branched evolution revealed by multiregion sequencing[END_REF][START_REF] Swanton | Intratumor heterogeneity: evolution through space and time[END_REF]. An alternative way to obtain the emergence of intra-tumor heterogeneity has been proposed in [START_REF] Lavi | The role of cell density and intratumoral heterogeneity in multidrug resistance[END_REF], by considering sufficiently large mutations in the models from [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF]. At this stage, let us stress that both structure variables, x and r, have a well defined biological meaning, so that they can be evaluated by means of laboratory experiments. In particular, a cell resistance level can be measured either by the expression level of ABC transporter genes, that are known to be associated with resistance to the drug [START_REF] Scotto | Transcriptional regulation of abc drug transporters[END_REF][START_REF] Szakács | Targeting multidrug resistance in cancer[END_REF], or by the minimal dose of each drug under consideration needed to kill a given percentage of the cell population [START_REF] Zhou | Sequential emergence of MRP-and MDR1-gene over-expression as well as MDR1-gene translocation in homoharringtonine-selected K562 human leukemia cell lines[END_REF].

Let us furthermore mention that the derivation of models able to include both selection and spatial dynamics, as the one presented here, is a key step toward a better comprehension of those mechanisms that underlie the evolution of ecological systems in general. These models can pave the way for interesting mathematical questions; see for instance [START_REF] Arnold | Existence of nontrivial steady states for populations structured with respect to space and a continuous trait[END_REF][START_REF] Bouin | Invasion fronts with variable motility: phenotype selection, spatial sorting and wave acceleration[END_REF][START_REF] Mirrahimi | Asymptotic analysis of a selection model with space[END_REF][START_REF] Mirrahimi | Dynamics of sexual populations structured by a space variable and a phenotypical trait[END_REF] and references therein.

The paper is organized as follows. In Section 2, we describe the mathematical model and the related underlying assumptions. Section 3 establishes the setup for numerical simulations and provides some considerations about the choice of the parameter functions and their biological consistency. Further details about the values of parameters can be found in Appendix B. Section 4 is devoted to study cell environmental adaptation in the framework of our model (i.e., how tumor cells adapt to the surrounding environment defined by nutrients and anti-cancer drugs). In particular, the evolution of phenotypic heterogeneity and chemotherapeutic resistance are analyzed in the presence of cytotoxic and cytostatic drugs. A qualitative mathematical justification for the results presented here is proposed in Appendix A. With the aim of supporting the design of optimal therapeutic strategies, in Section 5 we test the efficacy of therapeutic protocols based on constant infusion and/or bang-bang delivery (i.e., infusion schedules relying on bangbang control) of cytotoxic drugs, cytostatic drugs or combinations of cytotoxic and cytostatic drugs. Conclusions are drawn in Section 6, which provides also some ideas about future research perspectives.

A structured population model for a tumor cell spheroid exposed to anti-cancer drugs

We present a mathematical model for the dynamics of an in vitro population of tumor cells exposed to cytotoxic and cytostatic drugs. The cell population is assumed to be organized in a non-vascularized and radially symmetric spheroid, and to be structured by two non-negative real variables x ∈ [0, 1] and r ∈ [0, 1]. The former represents the normalized expression level of a cytotoxic resistant phenotype (i.e., roughly speaking, the level of resistance to cytotoxic agents), while the latter stands for the planar distance of cells from the center of the spheroid, whose radius is assumed to be normalized in order to have unitary length. Cancer cells proliferate through the consumption of nutrients -oxygen and glucose, since in this setting we do not consider the glycolytic phenotype (i.e., we do not distinguish between these two nutrients) -and die due to interpopulation competition for space and resources. We sum up the effects of cytotoxic drugs directly on mortality (i.e., in this simple setting, not involving the cell division cycle, we do not consider drug effects on cell cycle phase transitions [START_REF] Kimmel | Control theory approach to cancer chemotherapy: Benefiting from phase dependence and overcoming drug resistance[END_REF]), while we let cytostatic drugs slow down cellular proliferation. Nutrients and drugs are assumed to be in quasi-stationary equilibrium (i.e., their dynamics is assumed to be much faster than the selection dynamics involving cancer cells), and their concentrations vary over r due to diffusion, decay and consumption by cancer cells.

The population density of cancer cells is modeled by the function n(t, r, x) ≥ 0, so that local and global population densities at time t ∈ [0, T ] are computed, respectively, as

̺(t, r) = 1 0 n(t, r, x)dx, ̺ T (t) = 1 0 ̺(t, r)r 2 dr,
while the local average level of resistance X(t, r), the global average level of resistance χ(t) and the related variance σ 2 (t) can be evaluated as

X(t, r) = 1 0 x n(t, r, x) ̺ T (t) dx, χ(t) = 1 0 X(t, r)r 2 dr, σ 2 (t) = 1 0 1 0 x 2 n(t, r, x) ̺ T (t) r 2 drdx -χ(t) 2 .
In the present mathematical framework, the function σ 2 (t) provides a possible measure for intratumor phenotypic heterogeneity at time t. The function s(t, r) ≥ 0 identifies the concentration of nutrients available to cells. The densities of cytotoxic and cytostatic drugs are described, respectively, by c 1 (t, r) ≥ 0 and c 2 (t, r) ≥ 0. The dynamics of functions n, s, c 1 and c 2 is ruled by the following set of equations

∂ t n(t, r, x) = p(x) 1 + µ 2 c 2 (t, r) s(t, r) -d̺(t, r) -µ 1 (x)c 1 (t, r) n(t, r, x), ( 2.1)
-α s ∆s(t, r)

+ γ s + 1 0 p(x)n(t, r, x)dx s(t, r) = 0, (2.2)
-α c 1 ∆c 1 (t, r) + γ c 1 + 1 0 µ 1 (x)n(t, r, x)dx c 1 (t, r) = 0, (2.3)
-α c 2 ∆c 2 (t, r) + γ c 2 + µ 2 1 0 n(t, r, x)dx c 2 (t, r) = 0, (2.4)
with zero Neumann conditions at r = 0 coming from radial symmetry, and Dirichlet boundary conditions at r = 1

s(t, r = 1) = s 1 , ∂ r s(t, r = 0) = 0, c 1,2 (t, r = 1) = C 1,2 (t), ∂ r c 1,2 (t, r = 0) = 0, (2.5)
In the above equations, ∆ stands for the Laplacian in polar coordinates and:

• The function p(x) models the proliferation rate of cells expressing the resistance level x due to the consumption of resources. The factor

1 1 + µ 2 c 2 (t, r)
mimics the effects of cytostatic drugs. The parameter µ 2 models the average uptake rate and effect on proliferation of these drugs.

• The parameter d is the average death rate of cells due to the competition for limited space.

• The function µ 1 (x) denotes the death rate of cells due to the consumption of cytotoxic drugs.

• The parameters α s , α c 1 and α c 2 model, respectively, the diffusion constants of nutrients, cytotoxic drugs and cytostatic drugs.

• The parameters γ s , γ c 1 and γ c 2 represent the decay rate of nutrients, cytotoxic drugs and cytostatic drugs, respectively.

The total delivered dose of anti-cancer drugs is computed as

T 0 C 1 (t) + C 2 (t) dt.
The equation (2.1) can be recast in the equivalent form

∂ t n(t, r, x) = R x, ̺(t, r), c 1,2 (t, r), s(t, r) n(t, r, x),
in order to highlight the role played by the net growth rate of cancer cells, which is described by

R x, ̺(t, r), c 1,2 (t, r), s(t, r) := p(x) 1 + µ 2 c 2 (t, r) s(t, r) -d̺(t, r) -µ 1 (x)c 1 (t, r).
The following considerations and hypothesis are assumed to hold:

• With the aim of translating into mathematical terms the idea that expressing cytotoxic resistant phenotype implies resource reallocation (i.e., redistribution of energetic resources from proliferation-oriented tasks toward development and maintenance of drug resistance mechanisms), we assume p to be decreasing

p(•) > 0, p ′ (•) < 0. (2.6)
• The effects of resistance to cytotoxic therapies are modeled by assuming function µ 1 to be decreasing

µ 1 (•) > 0, µ ′ 1 (•) < 0. (2.7)
Further considerations and details about the definitions of p(x) and µ 1 (x) are provided in the following section.

Setup for numerical simulations and model parametrization

Numerical simulations are performed in Matlab making use of an implicit-explicit finite difference scheme combined with a shooting method with 200 × 200 points on the square [0, 1] × [0, 1]. For all simulations, the interval [0, T ] is selected as time domain, with time step dt = 0.1. We choose the boundary and initial conditions to be

s(t, r = 1) = s 1 > 0, c 1,2 (t, r = 1) = C 1,2 (t), (3.1
)

n(t = 0, •, x) = n 0 (•, x) := C 0 exp(-(x -0.5) 2 /ε). (3.2)
The parameter s 1 stands for the infusion rate of nutrients; the positive real functions C 1,2 (t) model the infusion rates of cytotoxic/cytostatic drugs and are defined, case by case, according to the situation under investigation in Section 4 and Section 5. Since we are dealing with a tumor spheroid of small size and micro-tumors are often derived from a single cell clone, we set C 0 ≪ 1 and ε ≪ 1. In particular, the choice (3.2) mimics a biological scenario where most of the cells are characterized by the same intermediate level of resistance to therapies at the beginning of observations (i.e., the cell population is almost monomorphic).

We use linear functions to model the proliferation rate p(x) and the death rate µ 1 (x), which is consistent with the method of measuring selection via multiple regression, that can be adapted from ecology [START_REF] Brodie | Correlational selection for color pattern and antipredator behavior in the garter snake thamnophis ordinoides[END_REF][START_REF] Grothey | Defining the role of panitumumab in colorectal cancer[END_REF]. In more detail, we set

p(x) := a 1 + a 2 (1 -x), µ 1 (x) := b 1 + b 2 (1 -x). (3.3)
The above definitions satisfy the generic properties set in (2.6)-(2.7). The parameters a 1,2 and b 1,2 are constrained by the following conditions: ) require the values of the average proliferation rate and the average mortality rate due to the effect of cytotoxic drugs to be equal to some values coming from experimental evidence [START_REF] Corbett | Tumor induction relationships in development of transplantable cancers of the colon in mice for chemotherapy assays, with a note on carcinogen structure[END_REF][START_REF] Grothey | Defining the role of panitumumab in colorectal cancer[END_REF]. We set a 1 = 0.1 [day -1 ] and b 1 = 0.1 [day -1 ], so that the orders of magnitude of these parameters are consistent with the existing literature [?, ?, ?]. The values of the other parameters are selected with exploratory aim.

b 1 ∝ a 1 , b 2 ∝ a 2 , (3.4) 1 0 a 1 + a 2 (1 -x)dx = 0.231 [day -1 ] ⇒ a 2 = 0.4620 -2a 1 [day -1 ], (3.5)
1 0 b 1 + b 2 (1 -x)dx = 0.8 [day -1 ] ⇒ b 2 = 1.6 -2b 1 [day -1 ]. ( 3

Study of cell environmental adaptation and phenotypic heterogeneity

For the model described in the previous section, we now study how tumor cells adapt to the surrounding environment defined by nutrients and anti-cancer drugs. In particular, we select the distribution presented in Fig. 1 as initial datum in order to mimic a biological scenario where most of the cells are characterized by the same intermediate level of resistance to therapies at the beginning of observations (i.e., the cell population is almost monomorphic at time t = 0).

F0b-eps-converted-to.pdf

Cell dynamics without therapies

We begin our study by analyzing the dynamics of cancer cells without therapies (i.e., C 1,2 (t) := 0). The results obtained are summarized by Fig. 2(a), which shows how, in the absence of therapeutic agents, cells characterized by lower resistance levels and thus, due to assumption (2.6), by stronger proliferative potentials, are selected. At each position r, a different trait X(T, r) is favored (i.e., for each value of r, n(T, r, x) concentrates in a different point X(T, r)).

In the framework of our model, this result can be justified by noting that the concentration of resources varies along the radius of the spheroid (i.e., s(T, r) attains different values at each r, see solid lines in Fig. 3), and different densities of available nutrients imply the selection of different levels of ability to get resources, thus providing the basis for the emergence of polymorphism within the tumor cell population considered. This is in agreement with experimental evidence. In fact, it has been shown that growing spheroids display a proliferation gradient where only the external cells do proliferate [START_REF] Laurent | Multicellular tumor spheroid models to explore cell cycle checkpoints in 3d[END_REF]. These data are confirmed by in situ experiments where only tumor cells closed to vessels do proliferate in the absence of drugs [START_REF] Minchinton | Drug penetration in solid tumours[END_REF].

Cell dynamics under infusion of cytotoxic or cytostatic drugs

At first, we consider the effects that constant infusions of cytotoxic drugs induce on cell dynamics, i.e., we perform simulations setting C 1 (t) := C 1 > 0 and C 2 (t) := 0. Fig. 2(b) highlights how cytotoxic drugs promote a selective sweep toward resistant phenotypes. A polymorphic scenario arises at the end of simulations also in this case; in fact, a different level of resistance X(T, r) is selected at any level within the spheroid (i.e., for each value of r, n(T, r, x) concentrates in a different point X(T, r)). By analogy with the case without therapies, this is due to the fact that the concentrations of nutrients and cytotoxic drugs vary along the radius of the spheroid (i.e., s(T, r) and c 1 (T, r) attain different values for different values of r, see solid and dashed lines in Fig. 3).

In order to study how cancer cells respond to the on-off switch of cytotoxic drug infusion, we perform simulations setting

C 1 (t) := C 1 > 0, if t ∈ [0, 20] 0, if t ∈ (20, 60], (4.1)
and keeping C 2 (t) := 0. The above definition mimics a biological scenario where cytotoxic drugs are delivered in the time interval [0, T /3] only. Fig. 4 highlights the selection of higher levels of resistance during the infusion of cytotoxic drugs, i.e., on the time interval [0, T /3], and higher level of proliferative potential in the absence of cytotoxic drugs, i.e., on the time interval (T /3, T ]. In fact, when the infusion of cytotoxic drugs is stopped, more proliferative, and thus less resistant, cancer clones are favored [START_REF] Gatenby | A change of strategy in the war on cancer[END_REF][START_REF] Gatenby | Adaptive therapy[END_REF]. The switch from the selection for resistance to the selection for proliferative potential occurs in a progressive and continuous way, rather than through jumps in the distribution over the traits (see the right panel in Fig. 4). Finally, we analyze the dynamics of cancer cells under the effects of constant infusions of cytostatic drugs, i.e., we run simulations setting C 1 (t) := 0 and C 2 (t) := C 2 > 0 for any t ∈ [0, T ]. Fig. 2(c) shows how the cell distribution at the end of simulations is still close to the initial one (see Fig. 1), that is, cytostatic drugs tend to slow down selection dynamics in the cancer cell population.

The results presented in this subsection lead us to conclude that phenotypic heterogeneity within solid tumors can come from cell adaptation to local conditions. In fact, cells characterized by different levels of proliferative potential and resistance to therapies can be selected depending on space position, in relation with the local concentrations of nutrients and anti-cancer drugs.

Considerations about intra-tumor heterogeneity

Here we study the evolution over time of the average level of resistance and the related variance. The results presented in Fig. 5 lead us to conclude that cytotoxic drugs cause an increase in time of the average level of resistance expressed by the whole cell population χ(t), while the related variance σ 2 (t) decreases. In the framework of our model, this result implies that cytotoxic drugs reduce intra-tumor heterogeneity. This is in agreement with the Gause competitive exclusion principle, and it is consistent with experimental observations for the fact that cytotoxic drugs reduce intra-tumor heterogeneity by favoring highly resistant cancer clones [START_REF] Ding | Clonal evolution in relapsed acute myeloid leukaemia revealed by whole-genome sequencing[END_REF][START_REF] Gottesman | Mechanisms of cancer drug resistance[END_REF][START_REF] Landau | Evolution and impact of subclonal mutations in chronic lymphocytic leukemia[END_REF][START_REF] Szakács | Targeting multidrug resistance in cancer[END_REF].

In conclusion to this study of cell environmental adaptation, we observe that cytotoxic drugs tend to limit diversity by favoring only resistant clones. This results in a bottleneck effect, which decreases intra-tumor heterogeneity.

Study of optimized therapeutic protocols

In this section, we compare the efficacy of different schedules of drug delivery with the aim of supporting the development of optimized therapeutic protocols. The effects of bang-bang infusion of cytotoxic or cytostatic drugs are compared with that of constant supply in Subsection 5.1, while the same kind of comparison for cytotoxic and cytostatic drugs in combination is provided in Subsection 5.2. Finally, the effects of therapeutic strategies that combine constant delivery of cytotoxic drugs with bang-bang infusion of cytostatic drugs, and vice-versa, are inspected in Subsection 5.3. The infusion schemes (i.e., the boundary conditions C 1,2 (t)) in use throughout this section are summarized by Fig. 6, where the parameters C a,b,c,d model the instantaneous delivered doses. The values of these parameters are selected in such a way that the total delivered dose remains always the same.

Constant vs bang-bang infusion of cytotoxic or cytostatic drugs only

At first, we study the efficacy of therapeutic protocols relying on bang-bang delivery of cytotoxic drugs only, and we compare the results obtained with the ones of constant infusion. We perform simulations with C 2 (t) := 0 and C 1 (t) defined as in the center panel of Fig. 6. The results obtained are compared with the outcomes of simulations performed with C 2 (t) := 0 and C 1 (t) defined as in the left panel of Fig. 6.

Constant infusion of cytotoxic drugs leads to a temporary reduction of the population density of cancer cells, which is then followed by a relapse caused by the emergence of resistance, see Fig. 8(a). On the other hand, the bang-bang infusion scheme with the same total dose leads to a population evolution that is qualitatively the same as that observed in the absence of therapies, apart from a temporary reduction of the population density of cancer cells during drug delivery (see Fig. 8(b) and compare it with Fig. 7).

Then, we develop the same kind of analysis for cytostatic drugs only. We perform simulations with C 1 (t) := 0 and C 2 (t) defined as in the center panel of Fig. 6. The results obtained are compared with the outcomes of simulations performed with C 1 (t) := 0 and C 2 (t) defined as in the left panel of Fig. 6. As we already know from Section 3, the constant infusion of cytostatic drugs tends to slow down the selection of high proliferative potentials by comparison with the case without drugs, see Fig. 9(a) and compare it with Fig. 7. On the other hand, the dynamics of cancer cells under bang-bang delivery of cytostatic drugs is qualitatively the same as that observed in the absence of therapies, apart from a temporary reduction of the population density of cancer cells during drug delivery (see Fig. 9(b) and compare it with Fig. 7).

As a conclusion to this section, we notice that the constant infusion of cytotoxic drugs leads to a temporary reduction of the population density of cancer cells, which is followed by a relapse caused by the emergence of resistance. As noted in the previous section, the constant infusion of cytostatic drugs slows down the dynamics of cancer cells. On the other hand, the bang-bang infusion of cytotoxic or cytostatic drugs weakly affects the dynamics of cancer cells by comparison with the case without therapies, apart from temporary reductions of the population density of cancer cells during drug delivery.

With the doses used in our tests, neither constant nor bang-bang infusion of only cytotoxic (cytostatic) drugs will push cancer cells toward extinction.

Constant vs bang-bang infusion of cytotoxic and cytostatic drugs

This subsection aims at making a comparison between the therapeutic effects of constant and bang-bang delivery of cytotoxic and cytostatic drugs in combination. Therefore, we perform simulations with C 1,2 (t) defined as in the center or the left panel of Fig. 6, respectively.

Bang-bang infusion leaves qualitatively unaltered the dynamics of cancer cells observed in the absence of therapeutic agents, apart from temporary reductions of the population density of cancer cells during drug delivery (see Fig. 10(b) and compare it with Fig. 7). On the other hand, the constant infusion scheme under consideration pushes cancer cells toward extinction, see Fig. 10(a). This is in agreement with the conclusions drawn in [START_REF] Lorz | Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies[END_REF] and consistent with experimental observations suggesting that combination therapies can be more effective [START_REF] Tabernero | Phase II trial of cetuximab in combination with fluorouracil, leucovorin, and oxaliplatin in the first-line treatment of metastatic colorectal cancer[END_REF][START_REF] Ye | Randomized controlled trial of cetuximab plus chemotherapy for patients with KRAS wild-type unresectable colorectal liver-limited metastases[END_REF]. Therefore we are led to conclude that, keeping equal the total delivered dose of drugs, if cytotoxic and cytostatic agents are used in combination, protocols relying on simultaneous bang-bang infusion can be less effective than protocols relying on simultaneous constant infusion.

As a conclusion to this section, we observe that effective anti-cancer treatments can be designed by making use of proper combinations between cytotoxic and cytostatic drugs. If these drugs are delivered together, constant supply is more effective than bang-bang infusion and will push cancer cells toward extinction.

Mixed constant/bang-bang infusions of cytotoxic and cytostatic drugs

Finally, using the two types of drugs in combination, we test the effects of delivery schedules relying on constant infusion of cytotoxic drugs and bang-bang infusion of cytostatic drugs, or vice-versa. Therefore, we perform simulations with C 1,2 (t) as in the right panel of Fig. 6. We define C 1 (t) as the fixed line and C 2 (t) as the dashed line, or vice-versa.

Bang-bang infusion of cytostatic drugs and constant infusion of cytotoxic drugs allows to control the tumor size, although a detectable number of cancer cells survives within the population, see Fig. 11(a). On the other hand, in good qualitative agreement with experimental observations [START_REF] Foo | Effects of pharmacokinetic processes and varied dosing schedules on the dynamics of acquired resistance to erlotinib in EGFRmutant lung cancer[END_REF], therapeutic protocols relying on bang-bang infusion of cytotoxic drugs and constant delivery of cytostatic drugs can push cancer cells toward extinction, see Fig. 11(b).

These results, together with that presented in Subsection 5.2, support the idea that more effective therapeutic protocols can be designed by using cytotoxic and cytostatic drugs in combination, with constant delivery for both drugs, or bang-bang infusion for cytotoxic drugs and constant infusion for cytostatic drugs. This can be justified by observing that cytostatic drugs in constant infusion slow down the emergence of resistance. As a consequence, the efficacy of cytotoxic drugs is enhanced, so that the same total dose of cytotoxic drugs can push cancer cells toward extinction.

As a conclusion to this section, we notice that therapeutic protocols relying on bang-bang infusion of cytotoxic drugs -constant delivery of cytostatic drugs are more effective than therapeutic protocols based on bang-bang infusion of cytostatic drugs -constant delivery of cytotoxic drugs.

The former make possible to push cancer cells toward extinction, while the latter allow only a control on tumor size.

Conclusions and perspectives

Departing from theories derived in other contexts of population biology and Darwinian-like selection, we have developed a space and phenotype structured population model for the dynamics of cancer cells exposed to cytotoxic and cytostatic drugs. Relying on the assumption that cells are organized in a radially symmetric spheroid, the present model takes explicitly into account the dynamics of resources and anti-cancer drugs, which define the cellular environment. In the present model, space structure together with diffusion of nutrients and therapeutic agents are the key ingredients providing the basis for the emergence of intra-tumor heterogeneity (i.e., the simultaneous selection of several levels of resistance/proliferative potential within the cancer cell population).

Study of cell environmental adaptation and phenotypic heterogeneity

In the framework of this model, we have first made use of numerical simulations to analyze the evolution of phenotypic heterogeneity and the emergence of resistance to therapies (see Section 3), and we have reached the following conclusions:

• Phenotypic heterogeneity within solid tumors can be explained, at least partially, by cell adaptation to local conditions. In fact, cells characterized by different levels of proliferative potential and resistance to therapies are selected depending on space position, in relation with the distributions of resources and anti-cancer drugs.

• Cytostatic drugs can slow down tumor evolution, while cytotoxic drugs can favor the selection of highly resistant cancer clones and cause a decrease in the heterogeneity with respect to the resistance trait. In the framework of our model, this is not in contradiction with the Gause competitive exclusion principle.

Study of optimized therapeutic protocols

As a second step, we have tested, in silico, the capability of different therapeutic protocols to effectively contrast cancer progression. The cases of constant infusion, bang-bang delivery and mixed constant -bang-bang infusion of the anti-cancer drugs under study have been considered (see Section 4). The results we have obtained lead us to the following conclusions:

• Looking for combination therapies relying on cytotoxic and cytostatic drugs is a more effective strategy for fighting cancer rather than using high doses of cytotoxic or cytostatic drugs only.

• Therapeutic protocols relying on constant supply of cytotoxic and cytostatic drugs, or bangbang infusion of cytotoxic drugs and constant delivery of cytostatic drugs are more effective than the other protocols considered here, since they can push cancer cells toward extinction. This can be justified by noting that cytostatic drugs slow down the growth of resistant cells and, thus, cytotoxic drugs can kill more cancer cells in the same time window.

• Since constant supplies of cytotoxic drugs might be excluded in practice for toxicity reasons, the bang-bang delivery of cytotoxic drugs associated with the constant infusion of cytostatic drugs might turn out to be optimal.

Perspectives

Future researches will be addressed to extend the present model in order to include the dynamics of healthy cells, and study how to handle unwanted toxic side effects related to the delivery of anti-cancer agents. In view of this, a possible research direction is to approach the toxicity problem at stake as an optimal control problem, where the systemic cost for human body related to cancer growth and therapy infusion has to be minimized under the constraint that a minimal concentration of healthy cells should be preserved inside the system. In the same way, in view of adaptive therapy [START_REF] Gatenby | Adaptive therapy[END_REF], preserving a minimal proportion of cancer cells under a sensitivity threshold x ≤ x (to be tuned) is another possible constraint to be represented in an optimal control problem. This general direction of research, using optimal control settings, clearly aims at therapeutic optimization in the clinic of cancers.

From a modeling point of view, the present theoretical study has taken into account (additionally to the space variable r) only one drug resistance trait for two different classes of anticancer drugs. This study should also be completed in the future by further work involving a multidimensional structure variable x including different resistant phenotypes to different classes of drugs. Furthermore, it might be interesting to consider other phenotypes related, for instance, to epithelial-to-mesenchymal transition, glycolytic metabolism and dormancy, which are involved in the emergence of cell subpopulations less accessible to treatment than the classical forms of in situ proliferating epithelial solid tumors.

From a more biological point of view, understanding what are the cell population characters (genetically or epigenetically determined?) aggregated in this structure variable (r, x) standing for spatial heterogeneity and quantitated drug resistance, and how it can be related in experimental measurements with gene expression or epigenetic enzyme activity, is a big challenge that must be faced in transdisciplinary studies gathering mathematicians and biologists. We know from Luria and Delbrück's princeps experiment [START_REF] Luria | Mutations of bacteria from virus sensitivity to virus resistance[END_REF] that some drug resistance, due to stochastic genetic mutations in cell populations, occurring prior to drug exposure, is likely to exist in cancer cell populations, all the more so as genome instability is a common feature of these cells [START_REF] Beketic-Oreskovic | Decreased mutation rate for cellular resistance to doxorubicin and suppression of MDR1 gene activation by the cyclosporin PSC 833[END_REF][START_REF] Inukai | Presence of epidermal growth factor receptor gene T790M mutation as a minor clone in non-small cell lung cancer[END_REF]. Finding out what are the respective parts played in drug resistance by purely stochastic processes [START_REF] Gupta | Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells[END_REF] on the one hand, and by more Lamarckian-like phenomena resulting from adaptation of the cells surviving a massive drug infusion, involving epigenetic mechanisms in response [START_REF] Sharma | A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations[END_REF], and thus justifying the use of deterministic models, is another challenge that we intend to tackle in forthcoming studies.

Appendix A. A qualitative mathematical justification for phenotypic selection

From a mathematical standpoint, using the considerations drawn in [START_REF] Lorz | Dirac mass dynamics in a multidimensional nonlocal parabolic equation[END_REF][START_REF] Mirrahimi | Asymptotic analysis of a selection model with space[END_REF][START_REF] Perthame | Dirac concentrations in Lotka-Volterra parabolic PDEs[END_REF], the long term dynamics of X(t, r) can be formally characterized by evaluating lim t→∞ R x = X(t, r), ̺(t, r), c 1,2 (t, r), s(t, r) = 0.

In order to verify it numerically, we write the following identity R X(T, r), ̺(T, r), c 1,2 (T, r), s(T, r) = 0, which implies, due to the definitions of the functions p and µ 1 provided in Appendix B,

X(T, r) = 1 + µ 2 c 2 (T, r) a 2 s(T, r) -b 2 c 1 (T, r)[1 + µ 2 c 2 (T, r)] s(T, r) 1 + µ 2 c 2 (T, r) (a 1 + a 2 ) -d̺(T, r) -(b 1 + b 2 )c 1 (T, r) .
Fig. 12 shows how the curves X(T, r) obtained from the formula above are in good agreement with the positions of the maximum points x M (T, r) of n(T, r, x)/̺(T, r) that result from numerical simulations.

Appendix B. Values of the parameters

. 6 )

 6 Relations(3.4) represent a possible way to model the relation between the rate of regression under chemotherapy and the rate of tumor growth established by the Norton-Simon hypothesis[START_REF] Simon | The Norton-Simon hypothesis: designing more effective and less toxic chemotherapeutic regimens[END_REF]. Conditions (3.5)-(3.6

  Subsection 4.1 deals with cell dynamics without drugs, while in Subsection 4.2 we analyze the effects of constant infusions of cytotoxic and cytostatic drugs. Considerations about the evolution of intra-tumor heterogeneity are drawn in Subsection 4.3. A qualitative mathematical justification for phenotypic selection is provided by Appendix A. Details about the values of parameters are provided in Appendix B.

Figure 1 :

 1 Figure 1: Initial cell distribution. Plots of

Figure 2 :

 2 Figure 2: T = 6000, α s = 0.08, α c 1 = 0.08, α c 2 = 0.2, γ s = 1, γ c 1 = 1, γ c 2 = 1, µ 2 = 800, d = 0.3620; (a) s 1 = 0.1, C 1 (t) := 0, C 2 (t) := 0; (b) s 1 = 0.3, C 1 (t) := 0.07, C 2 (t) := 0; (c) s 1 = 0.3, C 1 (t) := 0, C 2 (t) := 0.07.

Figure 3 :

 3 Figure 3: T = 6000, α s = 0.08, α c 1 = 0.08, α c 2 = 0.2, γ s = 1, γ c 1 = 1, γ c 2 = 1, s 1 = 0.3, µ 2 = 800, d = 0.3620; (left) C 1 := 0.07, C 2 (t) := 0; (right) C 1 := 0, C 2 (t) := 0.07.

Figure 4 :

 4 Figure 4: T = 3000, α s = 0.08,α c 1 = 0.08, α c 2 = 0.2, γ s = 1, γ c 1 = 1, γ c 2 = 1, s 1 = 0.3, µ 2 = 800, d = 0.3620, C 1 = 0.1, C 2 (t) := 0.

Figure 5 :

 5 Figure 5: T = 3500, α s = 0.08,α c 1 = 0.4, α c 2 = 0.4, γ s = 1, γ c 1 = 1, γ c 2 = 1, s 1 = 0.3, µ 2 = 800, d = 0.3620; (dashed lines) C 1 (t) := 0, C 2 (t) := 0; (dotted lines) C 1 (t) := 0, C 2(t) := 0.12; (solid lines) C 1 (t) := 0.12, C 2 (t) := 0.

Figure 6 :

 6 Figure 6: T = 6000, α s = 0.08,α c 1 = 0.08, α c 2 = 0.2, γ s = 1, γ c 1 = 1, γ c 2 = 1, µ 2 = 800, d = 0.3620; (left) s 1 = 0.1, C 1 (t) := 0, C 2 (t) := 0; (center) s 1 = 0.3, C 1 (t) := 0.07, C 2 (t) := 0; (right) s 1 = 0.3, C 1 (t) := 0, C 2 (t) := 0.07.

Figure 8 :Figures 9 - 10 :

 8910 Figure 8: T = 6000, α s = 0.08,α c 1 = 0.08, α c 2 = 0.2, γ s = 1, γ c 1 = 1, γ c 2 = 1, s 1 = 0.3, µ 2 = 800, d = 0.3620, C 1 (t) := 0, C 2 (t) := 0.

Figures 11 :

 11 Figures 11: T = 6000, α s = 0.08, α c 1 = 0.08, α c 2 = 0.2, γ s = 1, γ c 1 = 1, γ c 2 = 1, µ 2 = 800, d = 0.3620, s 1 = 0.3, C a = 0.05, C b = 1.275.

Figures 12 :Figure 2 : 1 0

 1221 Figures 12: T = 6000, α s = 0.08, α c 1 = 0.08, α c 2 = 0.2, γ s = 1, γ c 1 = 1, γ c 2 = 1, µ 2 = 800, d = 0.3620, s 1 = 0.3, C c = 0.05, C d = 1.275.

Figure 3 :

 3 Figure 3: Distributions of resources and drugs. Plots of s(t, r) (solid line) and c 1 (t, r) (dashed line) at t = 120, for C 1 (•) := C 1 > 0, C 2 (•) = 0 (left panel). Plots of s(t, r) (solid line) and c 2 (t, r) (dotted line) at t = 200, for C 1 (•) := 0, C 2 (•) := C 2 > 0 (left panel). The plot of s(t, r) at t = 50 in the case C 1.2 (•) := 0 is analogous to that presented here. The unit of time is days. Details of the parameters used are provided in Appendix B.

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: On-off switch of cytotoxic drug infusion. Plots of n(t, r, x)/̺(t, r) at t = 20 (left-top panel), n(t, r, x)/̺(t, r) at t = 60 (left-bottom panel) and

Figure 7 :Figure 9 :Figure 10 :

 7910 Figure 7: Cell dynamics without drugs. Plots of 1 0 n(t, r, x)r 2 dr (left panel) and ̺ T (t) (right panel). In agreement with the results presented in Section 3, cells characterized by a low expression level of resistance to cytotoxic therapies (i.e., strong proliferative potentials) are selected and intra-tumor heterogeneity is high. The unit of time is days. All values are normalized with respect to the initial global population density. To be compared with Figs. 8-11.

Acknowledgments. The research leading to this paper was (partially) funded by the french "ANR blanche" project Kibord: ANR-13-BS01-0004. T.L. was supported by the Fondation Sciences Mathématiques de Paris and by a public grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (reference: ANR-10-LABX-0098), and by the FIRB project -RBID08PP3J


E Guilmineau
email: emmanuel.guilmineau@ec-nantes.fr

R Duvigneau

J Labroquère

Optimization of jet parameters to control the flow on a ramp

Keywords: Flow control, Optimization, Ramp

This study deals with the use of optimization algorithms to determine efficient parameters of flow control devices. To improve the performance of systems characterized by detached flows and vortex shedding, the use of flow control devices such as oscillatory jets, are intensively studied. However, the determination of efficient control parameters is still a bottleneck for industrial problems. Therefore, we propose to couple a global optimization algorithm with an unsteady flow simulation to derive efficient flow control rules. We consider as test case an backward facing step with a slope of 25 • , including a synthetic jet actuator. The aim is to reduce the time-averaged recirculation length behind the step by optimizing the jet blowing/suction amplitude and frequency

Introduction

Flow-control technology has been widely researched to improve the aerodynamic performance of transportation systems such as aircrafts and cars. The common objectives in many cases are to prevent massive flow separations. Actuator devices, such as synthetic jets or vortex generators, have proved their ability to modify the flow dynamics. However, the determination of efficient flow control parameters, in term of location frequency, amplitude, . . . , is tedious and highly problem dependent [START_REF] Donovan | Active flow control appiled to an airfoil[END_REF][START_REF] Ekaterinaris | Active flow control of wing separated flow[END_REF]. To overcome this issue, the numerical simulation of controlled flows is often considered to determine optimal control parameters, and the use of an automated optimization procedure is more and more observed [START_REF] Duvigneau | Optimization of a synthetic jet actuator for aerodynamic stall control[END_REF][START_REF] Zymaris | Optimal location for suction or blowing jets using the continuous adjoint approach[END_REF]. Several studies have shown that the simulation of controlled flows is a difficult task, since result may be dependent on the turbulence close used [START_REF] Garnier | Evaluation of the unsteady RANS capabilities for separated flows control[END_REF][START_REF] Labroqure | Impact of turbulence closures and numerical errors for the optimization of flow control devices[END_REF] but also on numerical errors [START_REF] Labroqure | Impact of turbulence closures and numerical errors for the optimization of flow control devices[END_REF].

Therefore, this paper presents numerical results in the context of the optimization of flow control parameters for the flow on a ramp.

Flow solver

The ISIS-CFD flow solver, developed by "Ecole Centrale de Nantes" and CNRS, solves the incompressible unsteady Reynolds-Averaged Navier-Stokes equations. This solver is based on finite volume method to build a spatial discretization for the transport equations.

The incompressible unsteady Reynolds-Averaged Navier-Stokes equations can be written (using the generalized form of Gauss' theorem) as:

∂ ∂t V ρdV + S ρ( - → U - - → U d ) • - → n dS = 0 (1a) ∂ ∂t V ρU i dV + S ρU i (U j -U dj ) • n j dS = S (τ ij n j -pn i )dS ( 1b
)
where V is the domain of interest, or control volume, bounded by a closed surface S moving at a velocity -→ U d with a unit outward normal vector -→ n

where n j is the j th component.

-→ U and p are respectively the velocity and pressure fields. τ ij are the components of the Reynolds stress tensor.

All flow variables are stored at the geometric center of arbitrary shaped cells. Volume and surface integrals are evaluated with second-order accurate approximations. The face-based method is generalized to two-dimensional or three-dimensional unstructured meshes for which non-overlapping control volumes are bounded by an arbitrary number of constitutive faces, that means cells can be polyhedral. A centered scheme is used for the diffusion terms, whereas for the convective fluxes, the Gamma Differencing Scheme (GDS) [START_REF] Jasak | High resolution nvd differencing scheme for arbitrarily unstructured meshes[END_REF] is used for this study. Through a Normalized Variable Diagram (NVD) analysis [START_REF] Leonard | Simple high-accuracy resolution program for convective modelling of discontinuities[END_REF], this scheme enforces local monotonicity and convection boundedness criterion. For more details, see Queutey and Visonneau [START_REF] Queutey | An interface capturing method for freesurface hydrodynamic flows[END_REF].

The velocity field is obtained from the momentum conservation equations, and the pressure field is extracted from the mass conservation constraint, or continuity equation, transformed into a pressure-equation. The pressure equation is obtained in the spirit of Rhie and Chow [START_REF] Rhie | A numerical study of the turbulent flow past an isolated aerofoil with trailing edge separation[END_REF]. Momentum and pressure equations are solved in a segregated manner as in the SIMPLE coupling procedure [START_REF] Issa | Solution of the implicitly discretized fluid flow equations by operator-splitting[END_REF].

A second-order backward difference scheme is used for time discretization. All spatial terms appearing in equations (1a) and (1b) are treated in a fully implicit manner. In this paper, the geometry is fixed. Therefore, the velocity -→ U d in the equations (1a) and (1b) is null.

In the case of turbulent flows, additional transport equations for modeled variables are discretized and solved using the same principles. Various turbulence closures are implemented. The method features sophisticated turbulence models: apart from the classical two-equation k-ε and k-ω models, the anisotropic two-equation Explicit Algebraic Reynolds Stress Model (EARSM), as well as Reynolds Stress Transport Models, are available, see [START_REF] Deng | Comparison of explicit algebraic stress models and second-order turbulence closures for steady flow around ships[END_REF] and [START_REF] Duvigneau | On the role played by turbulence closures in hull ship optimization at model and full scale[END_REF]. Recently, a Detached Eddy Simulation (DES) approach has been introduced, see [START_REF] Guilmineau | Numerical simulation with a DES approach for automotive flows[END_REF].

We consider actuations based on synthetic jets. The actuation is implemented as boundary conditions. The velocity is imposed on the jet boundary and is defined as:

U = U j A(X)sin(ωt)d j (2)
with A(x) a unit profile function, U j the amplitude, ω the angular frequency with ω = 2πf j where f j is the frequency of the synthetic jet, and d j the direction of the jet. In this study, A(X) is a sine squared function and d j is perpendicular to the boundary.

Optimizer

A surrogate-based optimizer, included in the FAMOSA optimization toolbox developed at INRIA by Opale Project-Team, is used for this study. It is based on the construction of a kriging model to describe the variation of the objective function value, for example the recirculation length or the drag coefficient, with respect to control parameters. Kriging models, also called Gaussian process models, belong to response surface models, that allow to predict a function value f at a given point x, on the basis of a set of known function values

F N = {f 1 , f 2 , • • • , f N } at some points X N = {x 1 , x 2 , • • • , x N }, that
are stored in a database. F N is assumed to be one realization of a multivariate Gaussian process [START_REF] Williams | Prediction with Gaussian processes: From linear regression to linear prediction and beyond[END_REF] which has a joint Gaussian distribution

p (F N |C N , X N ) = exp -1 2 F T N C -1 N F N (2π) N det(C N ) (3)
for any collections of inputs X N . C N is the N × N covariance matrix, whose elements C mn give the correlation between the function values f m and f n obtained at points x m and x n . We assume that these values are correlated, since they correspond to underlying physical phenomena. This is expressed in terms of a correlation function

k, i.e., C mn (x m , x n ) = cov(f m , f n ) = k(x m ,
x n ; Θ) with Θ a set of hyper-parameters that are calibrated to maximize the likelihood of the observations.

It can be shown [START_REF] Gibbs | Efficient implementation of Gaussian processes[END_REF] that the probability density for the function value at a new point is

p(f N +1 |(X N , F N ), C N , x N +1 ) ∝ exp - (f N +1 -fN+1 ) 2 2σ 2 f N +1 (4)
where

fN+1 = k T N +1 C -1 N F N , σ2 f N +1 = κ -k T N +1 C -1 N k N +1 (5)
with

κ = k(x N +1 , x N +1 ; Θ) and k N +1 = [k(x 1 , x N +1 ; Θ), • • • , k(x N , x N +1 ; Θ)] T .
Thus, the probability density for the function value at the new point x N +1 is also Gaussian with mean fN+1 and standard deviation σf N +1 . Therefore, the most likely value at the new point x N +1 is fN+1 . This value will be considered as the prediction of the kriging model. The variance σ2

f N +1 can be interpreted as a measure of uncertainty in the value prediction. The function value can expected to vary in some range like [ f -3σ, f + 3σ].

The optimization strategy used in this study is based on the iterative construction of a kriging model. The use of such a model for optimization must be an iterative process since it is not possible to build a model accurate enough to find the optimal parameters in only one step. The model should be updated with the results from flow simulation until some convergence criterion is fulfilled. Therefore, the algorithm is organized in two phases. During the first one, an initial a priori database is constructed, that gathers the flow response (objective function values) corresponding to different control parameter values. The control parameters are chosen in order to explore uniformly the search space, according to a Design Of Experiments (DOE) method. During the second phase, a kriging model is constructed on the basis of available data and is used to determine which flow simulations should be carried out and added into the database. This second phase is then repeated until convergence of the algorithm, as illustrated in Figure 1.

The robustness and efficiency of this algorithm depends critically on the choice of the new simulations to be carried out during the second phase. To obtain satisfactory results, the capability of a kriging model to predict the uncertainty related to value predictions is used. Indeed, for a given design vector x, a kriging model provides not only a prediction of the function value f (x), but also an estimation of the uncertainty of this prediction σ(x). Therefore, one can use this knowledge to determine not only the point that minimizes the model, but also areas for which the model is uncertain, i.e. of poor quality. Evaluating such points promotes the minimization of the cost function as well as the improvement of the kriging model. Details on the enrichment procedure can be found in [START_REF] Duvigneau | Kriging-based optimization applied to flow control[END_REF]

Test-case

The ramp is a backward facing step with an slope of 25 • , and is selected to test the optimization method on a flow control problem. This test-case is a configuration chosen by the French research group "Flow Separation Control", entitled in French "GDR Contrôle des décollements". This configuration is represented in Figure 2.

The step height h is equal to 100 mm and the slope of the step is equal to 25 degrees. Two positions of the synthetic jet are investigated. The first is located just after the step. This position is called configuration A. The second is just before the step and is called configuration B. In both cases, the jet direction is normal to the wall. The width of the slot is d = 0.5 mm, i.e. 0.05h.

The inlet of the computational domain is located at 11.8h of the corner of the step and the outlet at 25h. The mesh is generated by Hexpress TM , an automatic unstructured mesh generator. This software generates meshes containing only hexahedrons. The mesh contains 73502 nodes and 71375 cells.

At the inlet, a uniform velocity U ref = 20 m/s, which corresponds to a Reynolds number of Re = 1.23×10 5 with respect to the height of the ramp, h, is imposed as well as the top of the computational domain. At the outlet, the pressure is imposed. The walls have a no-slip condition and the boundary layer is refined at the wall to fulfill the y + < 1 criterion, as shown Figure 4.

Results

For all the numerical simulations, the turbulence model used is the k-ω SST model of Menter.

Non controlled configuration

A configuration without actuation is firstly studied. In this case, the jet boundary condition is replaced by a wall.

Figure 5 presents a comparison of the velocity profile between the experimental data [START_REF] Thacker | Présentation des premiers résultats expérimentaux sur le montage du GDR[END_REF] and numerical results at X/h = -2.4. We observe a good agreement between the experimental data and the numerical result. This comparison confirms that the inlet is sufficient far to get the correct boundary thickness.

Figure 6 shows the streamlines obtained by the numerical simulation and in experiments. A massive separation exists and any vortex shedding appears. The recirculation is characterized by a length, called L r . The numerical simulation predicts L r /h = 5.333. The center of the recirculation is located at X C /h = 2.659 and Y C /h = 0.455.

Controlled configuration

A synthetic jet is used to reduce the recirculation length. The period of the synthetic jet, called T j , is described with 200 time steps and the numerical simulation is carried out with 20 jet periods. The averaged flow is obtained by using the last five periods.

The optimization parameters considered in this study are restricted to the amplitude and the frequency of the synthetic jet. To construct the a priori database used to explore the design space, the parameters of the synthetic jet are supposed to vary in the following ranges:

• 15 m/s ≤ U j ≤ 100 m/s for the amplitude,

• 50 Hz ≤ f j ≤ 600 Hz for the frequency.

For both jet configurations, an initial of 16 flows is generated using an optimized Latin Hypercube Sampling (LHS), as illustrated in Figure 7.

The time history of the optimization process is represented in Figure 8 for both configurations. The first step of optimization loop gives the jet configuration which leads to the minimum recirculation.

The meta-models for both configurations are presented in Figures 9 and10, for the configuration A and the configuration B respectively. For the configuration A, the best parameters of the synthetic jet are f j = 77.330 Hz and U j = 63.760 m/s. So, the Strouhal number is St = f j × h/U j ≃ 0.387. These parameters deal to the moment coefficient C µ = d U j /L r U ref = 0.9% and the reduced frequency F + = f j L r /U ref = 2.06. For the configuration B, the best parameters are f j = 62.086 Hz, U j = 88.193 m/s and the Strouhal number St ≃ 0.310. for this configuration, the non-dimensional parameters are C µ = 1.8% and F + = 1.66. For both configurations, the best parameters are obtained for a frequency of the synthetic jet close to the boundary of the design space. We also observe that the meta-model begins to oscillate. The possible sources of oscillations are the interpolating model that becomes too complex when the new points close together are added to the database and the difficulty to maximize the kriging model.

Figure 11 presents the streamlines of the averaged flow obtained with the best parameters of the synthetic jet. With these jet parameters, the recirculation length decreases. For the configuration A, the length is L r = 4.004h while for the configuration B, L r is equal to 4.061h. In both cases, a reduction of the time-averaged recirculation of about 24% is observed.

Figure 12 presents a comparison of the pressure coefficient for the non controlled case and the pressure obtained with the best parameters of the synthetic jet for each position of the actuator. The presence of the synthetic jet manifests by a significant increase in the suction peak and an extended domain of adverse pressure gradient.

Figures 13 and14 give some snapshots of the flow on one synthetic jet period for the configuration A and the configuration B, respectively. For both cases, the flow is obtained with the best parameters of the actuator. The flow evolution is similar for both jet configurations. Three recirculation bubbles are created during a time period of the actuator, and for each time step, the bubble position is similar between the two configurations.

Figures 15 and16 present a zoom close the actuator for the previous snapshots. When the blowing is activated, a bubble is created, see Figure 15(a) and 16(a). However, this bubble is more large with the configuration B. When the jet velocity is maximum, see Figure 15(b) and 16(b), the initial bubble goes down slides down the ramp. According to the position of the synthetic jet, the flow is different. With the configuration A, another bubble is present above the first bubble. It is not the case with the configuration B. With this configuration, a separation is present in front of the jet. Then, when the suction is activated, the flow is similar.

Conclusions

In this study, the kriging-based optimization algorithm of FAMOSA has been coupled with the unsteady Reynolds-averaged Navier-Stokes solver ISIS-CFD. This approach has been applied to optimize the amplitude and frequency of a synthetic jet in order to shorten the time-averaged recirculation length behind a ramp. Two positions of the actuator have been investigated: one just after the corner of the ramp and the second just before. Each configuration with its own jet parameters leads to the same reduction of the averaged bubble. The best parameters for the configuration with the actuator after the corner of ramp are C µ = 0.9% and F + = 2.06 while for the configuration with the actuator is before the corner the parameters are C µ = 1.8% and F + = 1.66. However, these two results give a similar overview of the flow with small differences close the actuator.

Figure 5 :

 5 Figure 5: Comparison of the velocity profile at X/h = -2.4 for the non controlled configuration

Figure 6 :

 6 Figure 6: Comparison of the streamlines for the non controlled configuration

Figure 7 :Figure 11 :Figure 12 : 23 (Figure 13 :

 711122313 Figure 7: Illustration of the initial sampling generating using LHS

(a) Configuration A (b) Configuration B Figure 3: Position of the synthetic jet

Figure 4: Computational domain

(a) Configuration A (b) Configuration B Figure 8: Optimization history

List of Figures


Thomas Degueule

Olivier Barais

Arnaud Blouin

Benoit Combemale

The K3 Model-Based Language Workbench

We introduce K3, a model-based language workbench that eases the engineering of domain-specific languages. K3 features state-of-the-art facilities that increase modularity and reusability of software language artifacts to decrease their development costs. Aspect-oriented and executable metamodeling are supported through the K3AL action language. K3SLE provides software language engineering facilities such as model polymorphism and language inheritance, supported by a modeloriented typing system. We present the main components of K3, their integration into Eclipse, and the main research questions they tackle. Finally, we present the plan of the tool demonstration.

I. INTRODUCTION

Models have been mainly used as descriptive artifacts allowing engineers to easily master systems complexity and ease their comprehension. Recently, however, the trend behind model-driven engineering (MDE) shifts this point of view towards active models, where models become the main artifacts from which software components (e.g. code, documentation, compilers) are automatically generated. Model-driven technologies are increasingly used in the industry [START_REF] Whittle | The state of practice in model-driven engineering[END_REF] notably for the definition of domain-specific languages (DSL) [START_REF] Fowler | Domain-specific languages[END_REF].

While DSLs are gaining popularity, both their definition and tooling (e.g. checkers, compilers) still require significant development efforts that must be balanced with their limited number of users (by definition). K3 is a language workbench that attempts to ease the definition of DSLs and their associated tooling using MDE techniques. It features state-of-the-art facilities to speed up the development of DSLs through increased modularity and reusability in the definition of language artifacts. K3 seamlessly integrates with the different tools defined around the Eclipse Modeling Framework (EMF) 1 . K3 is distributed as an Eclipse bundle providing standard development services such as editors and compilers 2 .

In the next sections, we review K3AL and K3SLE, that compose K3, and show their relations within the Eclipse runtime framework, as depicted in Figure I.

II. K3AL: ASPECT-ORIENTED AND EXECUTABLE METAMODELING

K3AL is a DSL built on top of the Xtend programming language [START_REF] Efftinge | Xbase: implementing domain-specific languages for java[END_REF]. Leveraging the open-class mechanism [START_REF] Clifton | MultiJava: Modular open classes and symmetric multiple dispatch for Java[END_REF], it allows to "re-open" existing meta-classes of a metamodel to insert new features such as attributes, references, or operations. This mechanism serves as the basis for the definition of aspects on metamodels elements, enabling aspect-oriented modeling that helps to model complex software artifacts composed of intertwined and cross-cutting concerns [START_REF] Jézéquel | Model driven design and aspect weaving[END_REF]. K3AL defines a set of active annotations that can be used in plain Xtend files to express, for instance, aspects, pre/post conditions, or invariants. K3AL uses Xbase [START_REF] Efftinge | Xbase: implementing domain-specific languages for java[END_REF] as the expression language for defining operations bodies. Listing 1 shows the definition of a set of aspects used to weave executability into a simple finite-state machine metamodel depicted in the left part of Figure II. An aspect is a plain Xtend class consisting of a set of attributes and methods. The @Aspect annotation is used to specify the target meta-class in the target metamodel into which the attributes and methods are woven. In our example, the aspects are used to add a current state to the state machine and to specify the operational semantics of the execution through the definition of a set of methods that effectively runs a state machine model. Once the weaving done, the new features are statically available anywhere in the system.

III. K3SLE: A META-LANGUAGE FOR MODEL TYPING

The metamodel of a DSL can be seen as the definition of a group of related types, that is, a set of constraints over admissible graphs of objects (models). With that definition, 3, a model m is a graph of interconnected objects respecting the type constraints M T defined by its metamodel M M . However, this definition differs in the current mainstream modeling workbenches. A metamodel M M is typically implemented with a set of classes in a given programming language (e.g. Java in EMF). Then, a model is concretely made of set of instances of these classes. Consequently, the XML serialization provided by EMF explicitly embeds a URI which refers to M M . Thus, in current model-driven technologies a model conforms to one and only one metamodel corresponding to the one used to create the model. As a matter of fact, many of the proposed reuse mechanisms in MDE depend on concrete metamodels [START_REF] Kusel | Reuse in model-to-model transformation languages: are we there yet? Software & Systems Modeling[END_REF]. So, while most of the tools would be reusable for a family of close DSLs, it is not possible in practice.

MM MT m MT' « subtypeOf » « typedBy » « implements » « conformsTo » Fig. 3.

Typing artifacts in MDE

The K3SLE meta-language aims at solving these problems by abstracting the overly restrictive conformance relation with a typing relation allowing to manipulate a model through different DSLs [START_REF] Degueule | Leveraging Family Polymorphism in MDE[END_REF]. The model-oriented type system of K3SLE leverages model typing [START_REF] Steel | On model typing[END_REF] and family polymorphism [START_REF] Ernst | Family polymorphism[END_REF] to allow any model m conforming to M M to be seen in a polymorphic way as typed by any super type M T ′ of M T . This model polymorphism mechanism enables the definition of generic model-manipulation tools defined over a model type M T . Thus, any model m can be supplied to these generic tools, providing that its metamodel M M implements a model type M T ′ that is a subtype of M T .

The main constructs of K3SLE are model types, metamodels, and transformations. Listing 2 shows a simple K3SLE program where the transformation execute is reused over two XMI models conforming to two different metamodels. Basic inheritance and implementation relations are also described. Aspects defined with K3AL can also be directly manipulated.

IV. TOOL DEMONSTRATION

The K3 demonstration we propose motivates the challenges it tackles and illustrates our approach and the different technical K3 components. Starting from the abstract syntax of a language defined using EMF, attendees will be introduced to aspect-oriented modeling through the use of Xtend and K3AL. K3SLE will then be used to develop representative transformations and tools around the language definition. The benefits of model polymorphism will be demonstrated through the definition of language extensions.

Fig. 1 .

 1 Fig. 1. K3 technology stack

Fig. 2 .

 2 Fig. 2. Weaving executability into a finite-state machine metamodel