Lionel Roques

François Hamel

Mathematical analysis of the optimal habitat configurations for species persistence

We study a reaction-diffusion model in a binary environment made of habitat and non-habitat regions. Environmental heterogeneity is expressed through the species intrinsic growth rate coefficient. It was known that, for a fixed habitat abundance, species survival depends on habitat arrangements. Our goal is to describe the spatial configurations of habitat that maximise the chances of survival. Through numerical computations, we find that they are of two main types -ball-shaped or stripe-shaped. We formally prove that these optimal shapes depend on the habitat abundance and on the amplitude of the growth rate coefficient. We deduce from these observations that the optimal shape of the habitat realises a compromise between reducing the detrimental habitat edge effects and taking advantage of the domain boundary effects. In the case of an infinite-periodic environment, we prove that the optimal habitat shapes can be deduced from those in the case of a bounded domain.

Introduction

In the recent literature, the loss of habitat connectivity, also known as environmental fragmentation, has been more and more often cited as a major cause of species endangerment. In the IUCN Red List of Endangered Species [START_REF]IUCN Red List of Threatened Species. A Global Species Assessment[END_REF], it is shown to be one of the most important threats on the endangered species.

The main causes of these habitat losses and fragmentation are related to human activities. For instance, roads create barriers for many animals, and carnivores are particularly vulnerable to the resulting habitat fragmentation (see e.g. Treves et al. [START_REF] Treves | Predicting Human-Carnivore Conflict: a Spatial Model Derived from 25 Years of Data on Wolf Predation on Livestock[END_REF]). Agricultural and forestry activities are the key drivers of habitat loss affecting birds [START_REF]IUCN Red List of Threatened Species. A Global Species Assessment[END_REF]. In freshwater, dam construction is a well-established cause of ecosystem fragmentation and species extinction (an example is given by Morita and Yamamoto [START_REF] Morita | Effects of habitat fragmentation by damming on the persistence of stream-dwelling charr populations[END_REF]). Habitat fragmentation can also be the result of climate changes. Indeed warming up can reduce the size of the habitat patches for some mountain-top species; it is also responsible for the shrinking of lakes and rivers, which may disconnect sub-populations.

Though it seems natural that habitat loss is a major threat to species persistence, the effects of habitat fragmentation per se, i.e., without habitat loss, 1 depend on the considered ecosystem. Nevertheless, most theoretical studies involving single-species models (see below for details and references) predict a negative impact of fragmentation on species persistence. Indeed, the larger the habitat patches are, the more likely the sub-populations will be important and therefore less susceptible to internal and external extinction factors. These theoretical results are supported experimentally for some animals and plants [START_REF] Pimm | On the risk of extinction[END_REF][START_REF] Matthies | Population size and the risk of local extinction: empirical evidence from rare plants[END_REF]. However, a recent analysis of several empirical studies, by Fahrig [START_REF] Fahrig | Effects of habitat fragmentation on biodiversity[END_REF], shows that fragmentation per se can also have positive effects on biodiversity, which is notably explained in her paper by interspecific interactions. The long-standing debate on "single large or several small" reserves (see Simberloff [START_REF] Simberloff | The contribution of population and community biology to conservation science[END_REF] for a discussion and further references on this subject) illustrates the complexity of the issue of the effects of habitat fragmentation.

In addition to the number and area of the habitat patches, their shape can also have an influence on species persistence. For some species, the habitat's edge has a detrimental effect, for instance because of increased predation, or increased risk of accident (see Chalfoun et al. [START_REF] Chalfoun | Nest predators and fragmentation: a review and meta-analysis[END_REF] and Lovejoy et al. [START_REF] Lovejoy | Edge and other effects of isolation on Amazon forest fragments[END_REF]). Stripe shaped habitat patches have large edge/area ratios, compared to round shaped habitat patches. Thus, it seems natural to claim, like Diamond [START_REF] Diamond | The island dilemma: Lessons of modern biogeographic studies for the design of natural reserves[END_REF], that species for which the edge has a detrimental effect will be more likely to survive in the latter case of round shaped habitat patches.

The important question of the effects of habitat fragmentation has stimulated the development of numerous population models taking account of environmental heterogeneity. In individual based models (IBM) each individual interacts with its environment, and simulations are run until a general behaviour of the population can be observed. Typically, with this approach, With and King [START_REF] With | Extinction thresholds for species in fractal landscapes[END_REF] showed that the probability of dispersal success increased with habitat connectivity. In some situations, this kind of models can be advantageous since their implementation does not require sophisticated mathematical tools while they can incorporate precise information on the species behaviour. However, they remain limited in terms of general theoretical analysis [START_REF] Kareiva | Analyzing insect movement as a correlated random walk[END_REF][START_REF] Marsh | The form and consequences of random walk movement models[END_REF][START_REF] Gross | Individual-based modeling[END_REF]. Conversely, diffusion models, although they sometimes make oversimplifying assumptions about movement of real organisms, provide an excellent tool for obtaining theoretical and qualitative results on populations dynamics, which can lead to a better understanding of some ecological processes. The idea of modelling population dynamics with these models from molecular physics has emerged at the beginning of the 20th century, with random walk theories of organisms, introduced by Pearson and Blackeman [START_REF] Pearson | Mathematical Contributions of the Theory of Evolution, XV, A Mathematical Theory of Random Migration[END_REF]. Then, Fisher [START_REF] Fisher | The wave of advance of advantageous genes[END_REF] and Kolmogorov et al. [START_REF] Kolmogorov | Étude de l'équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique[END_REF] used a reaction-diffusion equation as a model for population genetics. Several years later, Skellam [START_REF] Skellam | Random dispersal in theoretical populations[END_REF] used this type of model to study biological invasions, and he succeeded to propose quantitative explanations of observations, in particular of spreads of muskrats throughout Europe at the beginning of 20th Century. Since then, these models have been widely used to explain spatial propagation or spreading of biological species (bacteria, epidemiological agents, insects, fishes, plants, etc., see the books [START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF][START_REF] Murray | Mathematical biology[END_REF][START_REF] Okubo | Diffusion and Ecological Problems -Modern Perspectives[END_REF] for review).

The reaction-diffusion models that we study in this article correspond to a natural extension of the Fisher model to heterogeneous environments. This kind of model has been introduced by Shigesada et al. [START_REF] Shigesada | Traveling periodic waves in heterogeneous environments[END_REF] and studied recently [START_REF] Cantrell | Diffusive logistic equations with indefinite weights: population models in disrupted environments[END_REF][START_REF] Shigesada | Biological Invasions: Theory and Practice, Oxford Series in Ecology and Evolution[END_REF][START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF][START_REF] Berestycki | Analysis of the periodically fragmented environment model : II -Biological invasions and pulsating travelling fronts[END_REF], but it already appears that the heterogeneous character of the coefficients plays an essential role. Indeed, for different spatial configurations of the environment, the modelled species can tend to extinction or survive, disperse or not, depending on the arrangement of the habitat, as illustrated by the work of Roques and Stoica [START_REF] Roques | Species persistence decreases with habitat fragmentation: an analysis in periodic stochastic environments[END_REF]. These models have been considered either in bounded domains Ω or in an infinite domain. This last case, Ω = R N , is especially helpful when spreading phenomena are studied since it generates the existence of travelling waves [START_REF] Berestycki | Analysis of the periodically fragmented environment model : II -Biological invasions and pulsating travelling fronts[END_REF][START_REF] Aronson | Multidimensional nonlinear diffusions arising in population genetics[END_REF][START_REF] Berestycki | Front propagation in periodic excitable media[END_REF][START_REF] Kinezaki | Modeling biological invasion into periodically fragmented environments[END_REF].

Cantrell and Cosner [START_REF] Cantrell | Spatial Ecology via Reaction-Diffusion Equations[END_REF] have established, in the particular case of a 1dimensional bounded domain with reflecting boundary condition, that in a binary environment composed of habitat and non-habitat patches, it was better to concentrate all the habitat at a boundary of the domain. Shigesada and Kawasaki [START_REF] Shigesada | Biological Invasions: Theory and Practice, Oxford Series in Ecology and Evolution[END_REF] proved a comparable result in the case of a 1-dimensional infinite domain, composed of periodic regions of habitat and non-habitat. They found that it was better for species persistence to concentrate the habitat patches at the centre of the period cell. Berestycki et al. [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF] then generalised these results to the N-dimensional case, in both bounded and infinite-periodic situations. To do so, they gave a necessary and sufficient condition for species survival, based on the instability of the state 0 (where no individuals are present). This condition is equivalent to the negativity of the smallest eigenvalue λ 1 of the linearised elliptic operator around the 0 state (such a criterion was also used in [START_REF] Cantrell | Spatial Ecology via Reaction-Diffusion Equations[END_REF], in the bounded domain case). Using this criterion, they found a way of rearranging the habitat that decreases the value of λ 1 . Thus, given an initial spatial configuration of the habitat, the rearranged and more aggregated configuration always provides better chances of species survival for the Shigesada et al. model, while it corresponds to the same habitat's area. When the space dimension N is larger than 1, among these aggregated configurations, we do not know which one minimises λ 1 .

In the 1-dimensional case, the question of the optimal shape of the habitat is already solved. As a result of the works [START_REF] Shigesada | Biological Invasions: Theory and Practice, Oxford Series in Ecology and Evolution[END_REF][START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF][START_REF] Cantrell | Spatial Ecology via Reaction-Diffusion Equations[END_REF], this shape is simply an interval. However, even in the 2-dimensional case, as far as we know, nothing is known about this optimal shape (even its existence is not established), except that, if it exists, the optimal shape has to be stable by the above-mentioned rearrangement process, as a consequence of the results in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF]. Mathematically, the problem is the same as minimising the eigenvalue λ 1 with respect to a coefficient of the equation, corresponding to habitat configuration. Thus it is a part of the vast area of eigenvalue minimisation problems, among which the most celebrated is surely the Rayleigh problem [START_REF] Rayleigh | The Theory of Sound[END_REF] of minimising the lowest frequency of a membrane. It says that, in the simple homogeneous Fisher's case, in bounded domain with lethal boundary and fixed area A, the principal eigenvalue λ 1 of the Dirichlet Laplacian is minimal when Ω is a ball. Thus if the modelled species survives in a domain of area A, then it automatically survives in a ball-shaped domain of area A, while in other domains of area A, the species may go to extinction. This result is a consequence of the classical Rayleigh-Faber-Krahn inequality [START_REF] Faber | dass unter allen homogenen Membranen von gleicher Fläche und gleicher Spannung die kreisförmige den tiefsten Grundton gibt[END_REF][START_REF] Krahn | Über eine von Rayleigh formulierte Minimaleigenschaft des Kreises[END_REF][START_REF] Krahn | Über Minimaleigenschaft der Kugel in drei und mehr Dimensionen[END_REF][START_REF] Lumiste | A Centenary[END_REF].

In this paper, we address the question of finding an optimal habitat configuration in bounded environments with reflecting (also called "Neumann") boundary conditions and in infinite-periodic environments. We prove the existence of an optimal configuration. Then, in the 2-dimensional case, we approximate it thanks to new numerical computations, which exhibit two thresholds in the habitat abundance. Between these two thresholds, the optimal configuration has the shape of a "stripe". These numerical results lead to new questions: (i) is there a relation between the shapes of the optimal configurations in the bounded and periodic cases? (ii) does the optimal configuration depend on the equation coefficients amplitude between the habitat and non-habitat regions? (iii) in the particular situation of an environment composed of stripes, are large habitat stripes separated with large gaps better than thin stripes with thin gaps? We solve these questions analytically, in any space dimension; however, for the sake of clarity, most of our statements are presented in the 2-dimensional case.

2 Model formulation and preliminary mathematical results

Formulation of the model

The classical reaction-diffusion model in homogeneous environments of Fisher and Kolmogorov et al. corresponds to the following equation,

∂u ∂t = D∇ 2 u + u (µ -νu) in Ω ⊆ R N . (2.1)
The unknown u = u(t, x) corresponds to the population density at time t and position x = (x 1 , . . . , x N ). The number D > 0 measures the dispersion rate and ∇ 2 stands for the spatial dispersion operator

∇ 2 u = N i=1 ∂ 2 u ∂x 2 i
, N being the spatial dimension. The constant real numbers µ > 0 and ν respectively correspond to the intrinsic growth rate of the modelled species and to its susceptibility to crowding.

The extension of this model to heterogeneous environments that we study in this paper (Shigesada et al. model) can be written as follows,

∂u ∂t = D∇ 2 u + u (µ(x) -ν(x)u), in Ω ⊆ R N . (2.2)
We assume that µ(x) takes two values µ + or µ -, depending on x, with µ + > µ -. Regions where µ(x) = µ + correspond to "habitat" and are favourable zones, whereas regions where µ(x) = µ -are less favourable for the species ("nonhabitat" regions). When µ -< 0, it means that the birth rate is less than the mortality rate in the non-habitat region. In this case, as proved in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF], the species may tend to extinction. The function ν is bounded, ν ≥ 0 in Ω and ν is not equal to 0 almost everywhere, which corresponds to the classical negative density-dependence.

We consider two kinds of domains: the bounded case, where Ω is a piecewise C 1 bounded domain; and the periodic case, with Ω = R N and the coefficients µ(x) and ν(x) are periodic in the sense that there exist L 1 , . . . , L N > 0 such that µ(

x 1 + k 1 L 1 , . . . , x N + k N L N ) = µ(x 1 , . . . , x N ) and ν(x 1 + k 1 L 1 , . . . , x N + k N L N ) = ν(x 1 , . . . , x N
) for all integers k 1 , . . . , k N ∈ Z and for all (x 1 , . . . , x N ) in R N . We denote L = (L 1 , . . . , L N ) and we say that µ and ν are L-periodic. The period cell C is defined by

C := [0, L 1 | × • • • × [0, L N ].
(2.3)

In the bounded case, we assume that we have Neumann boundary conditions: ∂u ∂n = 0 a.e. on ∂Ω, where n is the outward unit normal to ∂Ω.

Note that the study in R N with periodic coefficients cannot reduce to a study on a torus since the solution u of (2.2) is not periodic as soon as the initial condition u(0, x) is not periodic.

In the periodic case, C is defined by (2.3). In the bounded case, we set C := Ω. Moreover, in both cases, we assume that there exist two Borel sets

C + , C -⊂ C, s.t.    µ(x) = µ + if x ∈ C + , µ(x) = µ -if x ∈ C -, C = C + ∪ C -, (2.4)
so that C + corresponds to the habitat regions and C -to the non-habitat regions.

Definition 2.1 Let µ be a function satisfying (2.4). We say that C + (µ) := C + is the habitat configuration associated to the growth rate function µ.

A criterion for species persistence

The necessary and sufficient condition of species survival that we use in this work is based on the sign of an eigenvalue λ 1 .

This number λ 1 is the principal eigenvalue of the linear operator -D∇ 2µ(x) and can be defined as the unique real number such that there exists a function Φ > 0 which satisfies

L 0 Φ := -D∇ 2 Φ -µ(x)Φ = λ 1 Φ in Ω, Φ > 0 in Ω, (2.5)
either with L-periodicity in the periodic case, or with the boundary condition ∂Φ ∂n = 0 a.e. on ∂Ω in the bounded case (Φ is also unique up to multiplication by positive numbers, see Amann [START_REF] Amann | Fixed point equations and nonlinear eigenvalue problems in ordered Banach spaces[END_REF] for further details). As we can see on formula (2.5), λ 1 does not depend on the coefficient ν. When necessary, we shall distinguish the principal eigenvalues in the periodic and bounded cases by denoting them λ 1,p and λ 1,b respectively. Furthermore, in order to emphasise the dependency of λ 1 with respect to µ, we shall sometimes write λ 1 (µ) for λ 1 .

Note that, given a function µ satisfying (2.4), the set C + is uniquely defined.

Reciprocally, the "patch function" µ is uniquely defined by C + , µ + and µ -, and, in the periodic case, the periodicity conditions. Thus, for fixed values of µ+ and µ -, we note

λ 1 [C + ] := λ 1 (µ).
Since the operator L 0 is self-adjoint, it is well-known that λ 1 satisfies the following formula:

λ 1 (µ) = min ϕ∈G\{0} C D|∇ϕ| 2 -µ(x)ϕ 2 C ϕ 2 , ( 2.6)
where G is defined by

G := H 1 per = ϕ ∈ H 1 loc (R N
) such that ϕ is L-periodic} in the periodic case and G := H 1 (C) in the bounded case.

Given a bounded continuous nonnegative initial datum u(0, x) = u 0 (x), the behaviour of the solution u(t, x) of (2.2) in the periodic case is given in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF], Theorem 2.6. Namely, if λ 1 ≥ 0, then u(t, x) → 0 uniformly in R N , thus the species tends to extinction, whereas if λ 1 < 0, u(t, x) → p(x) (uniformly on all compacts sets of R N ), where p(x) is the unique positive and bounded stationary solution of (2.2) i.e. p satisfies

D∇ 2 p(x) + p(x) (µ(x) -ν(x)p(x)) = 0, x ∈ R N .
(

The existence, uniqueness and periodicity of p are proved in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF], Theorems 2.1 and 2.4. In the bounded case, a similar result holds. Its proof, which is more classical, can be viewed as a straightforward adaptation of the results in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF] (see also Berestycki [38] for an analysis in the case of Dirichlet boundary conditions). Moreover, in the bounded case, the convergence to the stationary state, u(t, x) → p(x) as t → ∞, is uniform in Ω.

Mathematical formulation of the optimisation problem

The habitat proportion h ∈ [0, 1] being fixed, we define the family of habitat configurations with proportion h by:

S h := {Borel sets C + ⊂ C, such that |C + | = h|C|}, (2.8)
where |E| denotes the Lebesgue measure of a set E. This habitat proportion h and the values µ + and µ -of the growth rate in the habitat and nonhabitat regions being fixed, our goal is to describe the optimal habitat configurations for species survival, that is, the sets C + ∈ S h which minimise λ 1 [C + ]. For given values of µ + and µ -, we denote the optimal value of λ 1 [C + ] by λ 1 := min

C + ∈S h λ 1 [C + ],
and by S h,µ + ,µ -the set of optimal habitat configurations:

S h,µ + ,µ -:= {C + ∈ S h , such that λ 1 [C + ] = λ 1 }.
(2.9)

Remark 2.2 For a fixed h value, λ 1 depends on µ + and µ -. Moreover, as we will see in the next sections, the set S h,µ + ,µ -may also depend on µ + and µ - as soon as N ≥ 2.

The existence of such optimal configurations is given by:

Theorem 2.3 For all µ + > µ -in R, and for all h ∈ [0, 1], the set S h,µ + ,µ -is not empty.
Proof. This result is proved rigourously in Appendix A.

When needed, we may distinguish in the sequel the sets of optimal habitat configurations in the bounded and periodic cases, denoting them by S b h,µ + ,µ - and S p h,µ + ,µ -, respectively. (2.4). This result, which says that the optimal configuration is of "bang-bang" type, is proved in Nadin [39].

Remark 2.4 Actually, the minimum of λ 1 (µ) among the functions µ in {µ ∈ L ∞ (C), µ -≤ µ ≤ µ + a.e. and 1 |C| C µ = hµ + + (1 -h)µ -} is attained on a "patch function" satisfying

A numerical method for finding an optimal habitat configuration in the bounded case

We consider here the two-dimensional bounded case with

C = [0, L 1 ] × [0, L 2 ].
The domain is covered by n disjoint rectangles, and C + is the union of n + of these rectangles (see Fig. 1). Moreover, we set n = n r × n c where n r is the number of rows and n c is the number of columns in C. We have numerically computed the value of λ 1 for different configurations of C + and C -. Among the huge number of possible configurations (C n + n = n!/(n + !(n -n + )!)), the optimal ones can be sought in a smaller subset of admissible configurations. Indeed, consider the configuration C +, * nr,nc obtained with the following algorithm: 1. on each row, the rectangles that belong to C + are moved to the left; 2. then, on each column the rectangles that belong to C + are moved down. The old and new habitat configurations C + and C +, * nr,nc have the same area and it follows from Theorem 6.3 in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF], that λ 1 [C +, * ] ≤ λ 1 [C + ]. Thus, we only have to look for the optimal configuration among the habitat which are stable through this algorithm, that is, habitats concentrated in a corner of the domain. The number of configurations N conf corresponds to the number of ways of writing the integer n + as a sum of n c integers, each one being smaller than or equal to n r . It can be computed in a similar way as the "partition function" (see e.g. [START_REF] Goldberg | Combinatorial Analysis[END_REF]). For our computations, we fixed As the mesh is refined i.e. when n r and n c increase, the "discrete" optimal habitat configurations C +, * nr,nc , obtained by the above process, converge to an optimal configuration C +, * ∈ S h= n + n ,µ + ,µ -. Indeed, let us define the distance between two subsets A 1 and A 2 of C by the measure of the symmetric difference of

C = [0, 1] × [0, 1], µ + = 10, µ -= -1, n r = n c =
A 1 and A 2 , d set (A 1 , A 2 ) = |(A 1 ∪ A 2 ) -(A 1 ∩ A 2 )|, (2.10)
and the distance between a subset A ⊂ C and the family of subsets S h,µ + ,µ -by

d set (A, S h,µ + ,µ -) = inf K∈S h,µ + ,µ - d set (A, K), (2.11)
then, we have the following

Proposition 2.5 As n r ,n c → +∞, d set (C +, * nr,nc , S h,µ + ,µ -) → 0 as n r , n c → +∞.
Proof. This result is proved in Appendix B.

Deducing the periodic optimal configuration from the bounded case

With a numerical method comparable to that of §2.4, using Theorem 2.11 in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF], we could compute an optimal habitat configuration in the periodic case. However, the computations would be lengthy whereas we can directly deduce it from the bounded case. Indeed,

Theorem 2.6 For fixed numbers µ + , µ -∈ R, with µ + > µ -, and h ∈ [0, 1], let C +, * b be an optimal habitat configuration in S b h,µ + ,µ -in the bounded case. Assume that C +, * b is stable by monotone rearrangement with respect to each variable, with C = [0, L 1 | × [0, L 2 ]. Let C +, *
C +, * p := C +, * b ∪ R 1 [C +, * b ] ∪ R 2 [C +, * b ∪ R 1 [C +, * b ]],
where R i is the reflection with respect to the axis

{x i = L i }, for i = 1, 2. Then C +, * p is an optimal periodic configuration in S p h,µ + ,µ -on C = [0, 2L 1 |×[0, 2L 2 ].
Proof. In the bounded case, the existence of an optimal configuration, stable by monotone rearrangement with respect to each variable, follows from the proof of Theorem 2.3. The remaining part of the proof is given in Appendix C. We have proved a comparable result in any space dimension N . Thus, using our numerical computations in the bounded case, and from Proposition 2.5, we can derive an approximation of the optimal configuration in the periodic case.

Results

Both numerical and analytical results are presented in this section, in the bounded and periodic cases.

The bounded case

Using the numerical method presented in §2.4, we obtained that, when the proportion of habitat h = n + n is small, the optimal habitat regions have the shape of a quarter of disc (Fig. 2). When this proportion is close to 1, these optimal configurations look like the complement of a quarter of disc. When h is not too small nor too large, the optimal configurations are stripe-shaped, and are aggregated along a side of the boundary of the domain.

In the particular case of striped configurations, we have proved the following result:

Theorem 3.1 Let µ + > µ -, and h ∈ [0, 1]. Assume that L 1 = L 2 ,

and let S s

h be the subset of S h made of striped configurations. Then min

C + ∈S s h λ 1 [C + ] is
attained for the striped configuration concentrated along the smallest side of the domain.

Proof. See Appendix D. This result is illustrated on Fig. 3. When the habitat shape C + is fixed a priori as a ball, it is proved in [START_REF] Harrell | On the placement of an obstacle or a well so as to optimize the fundamental eigenvalue[END_REF] that the optimal position of the ball has to be tangent to the boundary of the domain, as depicted in Fig. 4. However, this configuration is never optimal among all habitats of fixed area when the shape of the habitat is not fixed a priori:

Proposition 3.2 Let µ + > µ -, and h ∈ [0, 1]. Assume that C + ∈ S h
is tangent to a boundary of the domain, in the sense that there exist a point x 0 ∈ ∂C and a ball In the bounded and periodic cases, the configuration (b), where the habitat C + is aggregated along the smallest side of the domain provides lower values of λ 1 than the configuration (a), independently of the habitat proportion h, and of µ + and µ -. Figure 4: A configuration which is not optimal in the bounded case. The habitats regions C + are depicted in black.

B ⊂ C + such that x 0 ∈ ∂B ∩ ∂(C\C + ) • . Then λ 1 [C + ] > λ 1 . Equivalently, this means that C + ∈ S h,µ + ,µ -. 0 1 1 (a) n + = 15 0 1 1 (b) n + = 45 0 1 (c) n + = 75 0 1 1 (d) n + = 80 0 1 1 (e) n + = 82 0 1 (f) n + = 83 0 1 1 (g) n + = 85 0 1 1 (h) n + = 90 0 1 (i) n + = 120 0 1 1 (j) n + = 150 0 1 1 (k) n + = 180 0 1 (l) n + = 210
Proof. This is a consequence of Steps 2 and 3 of Appendix E. See Remark 4.1.

The periodic case

As in the bounded case, when the proportion of habitat is small, the optimal configuration of C + is disc-shaped (Fig. 5). Then, for a larger proportion, it takes the form of a stripe, and when the proportion of habitat approaches 1, it looks like the complement of a disc. When the habitat region is stripeshaped, and when the period cell is not a square, two types of configurations exist. Similarly to Theorem 3.1, we formally prove in Appendix D that the configuration with wide stripes of habitat and wide gaps is always better than the configuration with narrow stripes and narrow gaps between the stripes.

These computations show that the shape of the optimal habitat configu-

0 2 2 (a) h = 0.2 2 0 2 (b) h = 0.4 2 2 0 (c) h = 0.6 2 2 0 (d) h = 0.8
Figure 5: Approximation of an optimal habitat configuration in the periodic case, with a period cell of size 2 × 2, in terms of the habitat abundance h = |C + |/|C|, with µ -= -1, µ + = 10. The black regions correspond to habitat and the white regions to non-habitat. The figures are assumed to extend in all direction indefinitely. ration is highly dependent on the proportion of habitat, even if two types of optimal shapes (disc and stripe) are numerically observed. This result rules out the attractive conjecture of the optimality of the striped configurations for all habitat abundance. Actually, we can prove formally the following result:

Theorem 3.3 Assume that N = 2. Let µ + > µ -be fixed.
(a) The disc configuration is not always optimal: assume that h = π(min{L1,L2}) 2

4|C|

, and that

C + is a disc in S h . Then λ 1 [C + ] > λ 1 .
(b) The stripe configuration is not always optimal: there exist

L 1 , L 2 > 0 and h ∈ [0, 1] such that for all stripe C + in S h , λ 1 [C + ] > λ 1 .
Proof. Part (a) follows from Appendix E, Steps 2 and 3, and Part (b) follows from Appendix F. We proved similar results in any space dimension N ≥ 2.

From Theorem 2.6, it follows that these results are also true in the bounded case, with, for Part (a), a quarter of disc instead of a disc.

3.3 Dependence of the optimal shape with respect to the difference between the growing rates µ + and µ - Let the proportion of habitat h be fixed. It then turns out that the optimal habitat configuration does not depend directly on the values of µ + and µ -, but only on the difference b := µ + -µ -. Indeed, as a consequence of formula (2.6), for all constant k ∈ R, we obtain

λ 1 (µ -k) = λ 1 (µ) + k. (3.12)
Thus, we may note S h,b the set of optimal habitat configurations, with b = µ + -µ -, instead of S h,µ + ,µ -.

In dimension 1, when C := [0, L 1 ], the optimal habitat shape is in fact totally independent of µ + and µ -. Indeed, in the periodic case, and up to a shift of the period cell, there is an unique configuration where C + is connected, which is therefore the optimal one (see [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF]). A similar result holds in the bounded case, where the optimal habitat configuration is an interval situated at either end of [0, L 1 ].

However, in dimension 2 or higher, given the proportion of habitat h, we prove that the optimal shape of the habitat does truly depend on the amplitude b of the growth rate µ:

Theorem 3.4 Assume that N ≥ 2. There exists h ∈ [0, 1] such that, for all C + ∈ S h , there exists b > 0 such that C + ∈ S h,b .
Proof. This is proved in Appendix E. The proof relies on the fact that (i) when the amplitude tends to infinity, the optimal habitat in the periodic case tends to aggregate in a ball; (ii) from Theorem 3.3, Part (a), the ball-shaped configuration of Fig. 6 is never optimal whatever the positive number b is.

This result of dependence with respect to amplitude is also true in the bounded case. Nevertheless, the transition between disc-shaped and stripe-Figure 6: A periodic configuration which is not optimal. The habitats regions are depicted in black.

shaped configurations as the proportion of habitat increases seems to hold for a wide range of values of b, as observed numerically in dimension 2. As the amplitude b increases, the ball configuration becomes better than the striped configuration for a wider range of habitat proportion (Fig. 7).

Discussion

We have numerically derived the optimal configurations of the habitat in the case of a 2-dimensional patchy environment, constituted of cells of habitat and non-habitat regions, distributed on a grid. As the size of the elementary cells decrease, these numerical configurations converge to an optimal habitat configuration for fixed values of habitat abundance and growth rates. This configuration is a global minimiser of λ 1 , contrary to most of the numerical methods for minimising eigenvalues under geometrical constraints (see e.g. [START_REF] Chanillo | Symmetry Breaking and Other Phenomena in the Optimization of Eigenvalues for Composite Membranes[END_REF][START_REF] Oudet | Numerical minimization of eigenmodes of a membrane with respect to the domain[END_REF]). The shape of these optimal configurations was not known before for dimensions higher than 1.

In rectangular bounded domains of R 2 , with reflecting boundaries, we have observed a transition effect from "quarter of disc" to "striped" and to "complement of a quarter of disc" optimal configurations as the proportion of habitat increases. The same phenomenon occurs in infinite periodic environments, with discs instead of quarter of discs. The transition occurs abruptly, and the striped configurations remain optimal for a substantial range of habitat proportion.

Actually, we have proved that the optimal configurations in the bounded and periodic cases can be deduced from each other. This means that, in a period cell, the role played by the proximity of the habitat patches lying in the other period cells is similar to the boundary effects of the bounded case.

From Theorem 6.3 in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF], and as recalled in the introduction of this paper, λ 1 decreases, and thus the chances of survival increase, when the habitat is rearranged in a monotone way, with respect to each space variable. Thus, the optimal configurations are always aggregated against the domain boundary. This underlines the favourability of the boundary effects. From an animal point of view, these effects can be interpreted as follows: the individuals situated in the habitat which try to cross the boundary are rejected into the habitat; therefore, the average time spent in the habitat region by each individual is larger when the habitat is closer to a boundary. From a vegetal viewpoint less descendants disperse in the unfavourable regions. Conversely, the edge of the habitat seems to have detrimental effects. This can again be illustrated with Theorem 6.3 in [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF]. Indeed, in the rearranged configurations the habitat has a small perimeter compared to situations where the habitat is more "fragmented". Hence, the optimal configuration is a compromise between lessening the detrimental habitat edge effects, which means that the perimeter of the habitat has to be minimal and leads to disc shapes, and taking advantage of the domain boundary (or the proximity of other patches in the periodic case), leading to striped shapes. However, both phenomena cannot be optimised simultaneously. When the habitat abundance is low, the optimal configuration looks like a quarter of disc in the bounded case, and looks like a disc in the periodic case; in these situation, the edge of the habitat is reduced, but the population does not take full advantage of the domain boundary. When the habitat abundance increases, the optimal configuration becomes stripe-shaped, and is concentrated along a side of the boundary of the domain. Thus, in that case, the positive influence of the boundary effects is important, compared to the detrimental habitat edge effects. However, when the domain (or period cell in the periodic case) is not a square, we prove formally that the best striped habitat configuration is aggregated along the smallest side of the domain. Thus this configuration reduces the edge in contact with non-habitat regions, while it does not take full advantage of the boundary, still realising a compromise between these edge and boundary effects. In the periodic case, this means that the configuration with wide stripes is always better than the one with narrow stripes. This gives an element of explanation of the observation of Kinezaki et al. [START_REF] Kinezaki | Modeling biological invasion into periodically fragmented environments[END_REF] that, in a periodic environment composed of stripes, when the scale of fragmentation is enlarged without changing the habitat proportion, the spreading speed of the modelled species increases.

We have proved that this equilibrium between boundary and edge effects also depends on amplitude of the difference between the growth rates in the habitat and non-habitat regions. As this amplitude grows, the detrimental effects of the habitat edge become more important compared to the favourable boundary effects, resulting in ball-shaped optimal configurations in highly contrasted environments. Nevertheless, for a wide range of values of the growth rate function, we observed the same transition behaviour between disc-shaped and striped configurations as the habitat abundance increases.

Our results in the bounded case largely depend on the reflecting (Neumann) boundary conditions. These boundary conditions mean that individuals encountering the boundary are reflected inside the domain. They can be encountered in many real-world situations, at different scales, corresponding, for instance, to the walls of a laboratory arena for small terrestrial species (see e.g. the experiment of Hannunen and Ekbom [START_REF] Hannunen | Host plant influence on movement patterns and the subsequent distribution of Lygus rugulipennis (Heteroptera: Miridae) nymphs[END_REF]), to fences or rivers surrounding gardens, forests or natural parks [START_REF] Shigesada | Biological Invasions: Theory and Practice, Oxford Series in Ecology and Evolution[END_REF][START_REF] Jaeger | Effects of Road Fencing on Population Persistence[END_REF] for larger species. This type of boundary conditions can also correspond to coasts, in the case of aquatic species. They can also be used for modelling plant dispersal, particularly when the seeds are dispersed by animals. Of course, other boundary conditions could have been envisaged. Of particular practical interest are also the absorbing (Dirichlet) conditions, that we briefly discussed in the introduction section. These conditions, where the domain's boundary "absorbs" all individuals encountering it, are often envisaged for species with passive movement, such as plants. Although most of the methods developed in this paper still work under these conditions, we have not considered here the issue of finding an optimal habitat configuration in this case. Firstly, because some numerical and theoretical results already exist in the mathematical literature [START_REF] Chanillo | Symmetry Breaking and Other Phenomena in the Optimization of Eigenvalues for Composite Membranes[END_REF]. Moreover, the question of optimal configuration should be of less interest in this case. Indeed, due to the negative effect of the domain's boundary, the habitat will concentrate at the centre of the domain. Therefore we cannot expect great changes in the topology of the optimal habitat configuration with respect to the model parameters.

A natural question arises: what is the interest of minimising λ 1 as soon as λ 1 < 0, since species survival is already granted? The answer is partially given by (3.12): λ 1 (µ -k) = λ 1 (µ) + k for all k ∈ R. Thus, a species in an environment with smaller values of λ 1 will be more robust to perturbations in its habitat. Moreover, as shown by [START_REF] Chekroun | Models of population dynamics under the influence of external perturbations: mathematical results[END_REF] it is also more robust to external perturbations corresponding for instance to harvesting, in the sense that δ can be chosen higher in ∂u ∂t = D∇ 2 u + u (µ(x) -ν(x)u) -δ, without risking the species extinction. Another natural question is: when the environmental heterogeneity also concerns the motility of the individuals, what results can be expected? Consider the equation

u t = ∇ 2 (D(x)u) + u (µ(x) -ν(x)u), (4.13)
where the diffusion term ∇ 2 (D(x)u) is in the Fokker-Planck form. The coefficient D(x) corresponds to the species motility (see e.g. [START_REF] Turchin | Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants[END_REF][START_REF] Okubo | Diffusion and Ecological Problems -Modern Perspectives[END_REF]). Setting v(t, x) := D(x)u(t, x), and using some results of [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF], we can obtain that the equation (4.13) exhibits the same behaviour as (2.2), replacing the initial growth term µ(x) by µ(x)/D(x). Thus, our above optimisation results also hold for equation (4.13), when D(x) is rearranged along with µ(x), and in the same way. Such a rearrangement is biologically meaningful, since the motility and the growth terms are generally associated with a certain type of habitat, and are therefore spatially correlated. Note that, for most species that engage in active habitat selection, such as birds, mammals of amphibians, D(x) is all the smaller as µ(x) is large, since they tend to slow down as they encounter favourable areas (e.g. [START_REF] Hannunen | Host plant influence on movement patterns and the subsequent distribution of Lygus rugulipennis (Heteroptera: Miridae) nymphs[END_REF]). Thus the functions µ and µ/D vary in the same way.

For the sake of simplicity of the following calculations, the cell C has been centred at the origin. Let P h be the set of measurable functions defined by:

P h := {µ satisfying (2.4) with |C + | = h|C|}, (4.14)
and

λ 1 := inf µ∈P h λ 1 (µ) = inf C + ∈S h λ 1 [C + ],
where λ 1 (µ) denotes the principal eigenvalue of the operator -D∇ 2 -µ with L-periodicity in the variables x 1 , . . . , x N (each function µ in P h can be extended by L-periodicity). By definition of λ 1 , there exists a sequence

(µ n ) n∈N in P h , such that λ 1 (µ n ) → λ 1 as n → +∞. Since, for all n ∈ N, λ 1 (µ * n ) ≤ λ 1 (µ n ) and µ * n ∈ P h , we also have λ 1 (µ * n ) → λ 1 as n → +∞, where µ *
n is the symmetric decreasing Steiner rearrangement of µ n with respect to each variable x 1 , ..., x N in the cell C (see [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF] for details). Up to a change on a set of zero measure, the functions µ * n can be constructed so that, for each 1

≤ i ≤ N and (x 1 , . . . , x i-1 , x i+1 , . . . , x N ) ∈ [-L 1 /2, L 1 /2] × • • • × [-L i-1 /2, L i-1 /2] × [-L i+1 /2, L i+1 /2] × • • • × [-L N /2, L N /2], there exists a ∈ [0, L i /2] such that the function s → µ * n (x 1 , . . . , x i-1 , s, x i+1 , . . . , x N ) is equal to µ + in [-a, a] and to µ -in [-L i /2, L i /2]\[-a, a].
Up to the extraction of some subsequence, Define, for all ϕ ∈ H 1 per \{0} and µ ∈ P h ,

µ * n (x) → µ * (x) ∈ {µ -, µ + } as n → +∞, for all x in the countable set Q N ∩ C. Let Ω + be the closure of {x ∈ Q N ∩ C, µ * (x) = µ + }. We extend the function µ * to C by setting µ * (x) = µ + if x ∈ Ω + and µ * (x) = µ -otherwise; next,
Q µ (ϕ) = C D|∇ϕ| 2 -µ(x)ϕ 2 C ϕ 2 . (4.16)
Then,

Q µ * (ϕ) ≤ Q µ * n (ϕ) + C |µ * -µ * n |ϕ 2 C ϕ 2
for all ϕ ∈ H 1 per \{0}.

Call Φ n the principal eigenfunction associated to µ * n , with the normalisation condition ||Φ n || ∞ = 1. That is, Φ n is L-periodic, positive, and it satisfies

-D∇ 2 Φ n -µ * n Φ n = λ 1 (µ * n )Φ n in R N . (4.17)
By multiplying this equation by Φ n and integrating over C, we get

Q µ * n (Φ n ) = λ 1 (µ * n ). Hence, Q µ * (Φ n ) ≤ λ 1 (µ * n ) + C |µ * -µ * n | C Φ 2 n , and since λ 1 (µ * ) = min ϕ∈H 1 per \{0} Q µ * (ϕ), λ 1 (µ * ) ≤ λ 1 (µ * n ) + C |µ * -µ * n | C Φ 2 n . (4.18)
Moreover, since µ -≤ µ * n ≤ µ + in R N and -µ + ≤ λ 1 (µ * n ) ≤ -µ -, the normalisation condition Φ n ∞ = 1 together with standard elliptic estimates (see Gilbarg and Trudinger [48]) imply that the sequence (Φ n ) n∈N is bounded in C 1 (R N ) (at least), whence inf n∈N Φ n L 2 (C) > 0. Thus, passing to the limit as n → +∞ in (4.18), we get, since µ * n → µ * in L 1 (C), λ 1 (µ * ) ≤ λ 1 . It then follows from the definition of λ 1 that λ 1 (µ * ) = λ 1 . Since µ * ∈ P h , the associated habitat configuration C +, * , defined by (2.4), belongs to S h , and λ 1 [C +, * ] = λ 1 , therefore C +, * is an optimal configuration. The infimum in the definition of λ 1 is then a minimum.

B. Convergence of the computed configurations to optimal configurations as the mesh is refined

This time, the proof is done in the bounded case of a rectangle. We assume that From the result of Appendix A, there exists a function µ * ∈ P h such that λ 1 (µ * ) = λ 1 . Let Φ * ∈ H 1 (C) ⊂ L 2p (C) for some p > 1, be the eigenfunction associated to the eigenvalue λ 1 , with Neumann boundary conditions. Let P nr,nc h be the set of all elements σ nr,nc of P h such that C + (σ nr,nc ) belongs to the grid of size n = n r × n c . For all ε > 0, and for n r and n c large enough, it is easily seen that there exists µ nr,nc ∈ P nr,nc h such that ||µ * -µ nr,nc || L p ′ (C) < ε (with

C := [0, L 1 ] × • • • × [0, L N ] . Let µ + , µ -∈ R,
1 p + 1 p ′ = 1
). Let Φ nr,nc be the principal eigenfunction associated to µ nr,nc and λ 1 (µ nr,nc ). Using formula (2.6), and since, when µ = µ * , the minimum in this formula is precisely attained for Φ = Φ * , we obtain

λ 1 (µ nr,nc ) ≤ λ 1 (µ * ) + C |µ nr,nc -µ * |(Φ * ) 2 C (Φ * ) 2 ≤ λ 1 (µ * ) + ε Φ * 2 L 2p (C) Φ * -2 L 2 (C) , (4.19
) from Hölder inequality. Setting α

:= Φ * 2 L 2p (C) Φ * -2 L 2 (C) we get λ 1 (µ nr,nc ) ≤ λ 1 (µ * )+αε. Since P nr,nc h ⊂ P h , we therefore have λ 1 (µ nr,nc ) ∈ [λ 1 (µ * ), λ 1 (µ * )+ αε]. Let µ *
nr,nc be a function which minimises λ 1 in P nr,nc h ; up to a rearrangement we can assume that µ * nr,nc is nonincreasing with respect to each variable

x 1 , • • • , x N . Then λ 1 (µ * ) ≤ λ 1 (µ * nr,nc ) ≤ λ 1 (µ nr,nc ); thus, λ 1 (µ * nr,nc ) ∈ [λ 1 (µ * ), λ 1 (µ * ) + αε], and λ 1 (µ * nr,nc ) → λ 1 (µ * ) = λ 1 as n → ∞.
From this convergence result, and using the same arguments as in Appendix A, we know that, up to the extraction of some subsequence, µ * nr,nc → µ * ∞ in L 1 (C) as n r , n c → +∞, where µ * ∞ ∈ P h satisfies λ 1 (µ * ∞ ) = λ 1 . Let us define the set of "optimal functions" µ by:

P h := {µ ∈ P h , such that λ 1 (µ) = λ 1 }. (4.20)
Since we can construct a converging subsequence from every subsequence of (µ * nr,nc ) nr,nc∈N , the computed configuration µ * nr,nc is as close as we want to the set P h for n r , n c large enough, in the sense that

dist L 1 (C) (µ * nr,nc , P h ) → 0 as n r , n c → +∞. (4.21)
Let C +, * nr,nc be the habitat configuration associated to µ * nr,nc by (2.4). Then

dist L 1 (C) (µ * nr,nc , P h ) = |µ + -µ -|d set (C +, * nr,nc , S h,µ + ,µ -),
where S h,µ + ,µ -is defined by (2.9) and d set (C +, * nr,nc , S h,µ + ,µ -) is defined by (2.10-2.11). Since

µ + = µ -, d set (C +, * nr,nc , S h,µ + ,µ -) → 0 as n r , n c → +∞.
C. Deducing periodic optimal configurations from bounded optimal configurations From Appendix A, we know that there exists a function µ * b in P h , nonincreasing with respect to each variable

Let C = [0, L 1 ]ו • •×[0, L N ] , µ + , µ -∈ R, and h ∈ [0, 1]. Set C := [0, 2L 1 ]ו • •× [0, 2L N ]. Let
x 1 , • • • , x N . Let us show that µ * p (x 1 , • • • , x N ) := µ * b (|x 1 -L 1 |, • • • , |x N -L N |) minimises λ 1,
P h := {µ ∈ Ph , such that λ 1,p (µ) = λ 1,p }.
Let µ p be a symmetric nonincreasing function in P h and Ψ p be the corresponding symmetric nonincreasing eigenfunction (see [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF]). Set

Ψ b (x 1 , • • • , x N ) := Ψ p (x 1 + L 1 , • • • , x N + L N ).
The function Ψ b is well defined on C, and by symmetry and periodicity properties of Ψ p , it is of class C 1 (C) and it satisfies Neumann boundary conditions. Moreover, Ψ b > 0,

-D∇ 2 Ψ b -µ p (x 1 + L 1 , • • • , x N + L N )Ψ b = λ 1,p Ψ b in C,
and µ p (• + L 1 , • • • , • + L N ) ∈ P h . By definition of λ 1,b , we get λ 1,p ≥ λ 1,b .
Let Φ b be the eigenfunction corresponding to λ 1,b with the configuration µ * b . Let us define

Φ p (x 1 , • • • , x N ) := Φ b (|x 1 -L 1 |, • • • , |x N -L N |) in C, and let µ * p (x 1 , • • • , x N ) = µ * b (|x 1 -L 1 |, • • • , |x N -L N |).
Then, from the Neumann boundary conditions satisfied by Φ b , it follows that Φ p is a weak solution of

-D∇ 2 Φ p -µ * p (x)Φ p = λ 1,b Φ p in C with periodic boundary conditions. Since | C+ (µ * p )| = 1 µ + C∩{µ * p =µ + } µ * p = 2 N µ + C∩{µ * b =µ + } µ * b = 2 N |C + (µ * b )| = h| C|,
it follows that µ * p ∈ Ph . Thus, we get λ 1,p ≤ λ 1,b , and finally λ 1,p = λ 1,b . As a conclusion, we obtain that µ * p ∈ P h , and therefore that C+ (µ * p ) is an optimal habitat configuration in the periodic case, with period cell C. Note that C+ (µ * p ) simply corresponds to successive reflections with respect to each axis {x i = L i }, for i = 1 . . . N starting from the optimal configuration C + (µ * b ) in the bounded case.

D. Proof of the comparison principle between the striped configurations

The proof is done in the periodic case. Let λ V be the value of λ 1 in the case of the vertical configuration depicted in Fig. 3 (a), and λ H be the value of λ 1 for the horizontal configuration depicted in Fig. 3 (b). Then, a simple change of variable leads to λ

V [L 1 , L 2 ] = λ H [L 2 , L 1 ] (in λ V [X, Y ] and λ H [X, Y ],
X corresponds to the length of the period cell in the x 1 direction and Y corresponds to the length of the period cell in the x 2 direction). Since

L 1 L 2 = |C|, we get λ V [L 1 , |C|/L 1 ] = λ H [|C|/L 1 , L 1 ]. (4.22)
Let us now prove that λ V [L 1 , |C|/L 1 ] decreases with L 1 . Assume that we are in the configuration of the vertical stripe. Then µ(x 1 , x 2 ) does not depend on x 2 . Set μ(x 1 ) := µ(x 1 , 0). The function μ is L 1 -periodic and, in the period cell [0, L 1 ], μ is equal to µ + on an interval of length l 1 = |C + |/L 2 , and it is equal to µ -on an interval of length L 1 -l 1 . Let Φ and λ V = λ 1 be defined by equation (2.5). Actually, since Φ is uniquely determined by (2.5) and the periodicity conditions, it depends on x 1 only. Set ψ(x 1 ) := Φ(x 1 , 0). Then ψ and λ V satisfy the following one-dimensional equation: Then, since the operator ϕ → -(L 1 ) -2 Dϕ ′′ -μϕ is self-adjoint, we have the following formula:

-Dψ ′′ (x 1 ) -μ(x 1 )ψ(x 1 ) = λ V ψ(x 1 ), x 1 ∈ R, ψ > 0 in R (4.
λ V = min ϕ∈H 1 per \{0}     1 L 2 1 × 1 0 D(ϕ ′ ) 2 (y)dy 1 0 ϕ 2 (y)dy - 1 0 μ(y)ϕ 2 (y)dy 1 0 ϕ 2 (y)dy     , (4.25)
where H 1 per denotes here the set of H 1 loc (R) functions which are 1-periodic. From (4.25), we get that λ V = λ V [L 1 , |C|/L 1 ] is decreasing with respect to L 1 . Notice that the strict monotonicity follows from the fact that the minimum in (4. [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF]) is reached at the function ξ, which is not constant since μ is not constant. As a conclusion, using (4.22), it follows that λ

V [L 1 , L 2 ] > λ H [L 1 , L 2 ] if L 1 < L 2 .
E. Dependence of the optimal configuration on the amplitude of the variations of the environment

We consider here the periodic case. We fix N ≥ 2. In the sequel, for any Borel set C + ⊂ C • := (0, L 1 ) × • • • × (0, L N ) and for any real numbers µ -< µ + , we call

λ 1 [C + , µ + , µ -] = min ϕ∈H 1 per \{0} C D|∇ϕ| 2 - C µ(x)ϕ 2 C ϕ 2 (4.26)
the principal eigenvalue of the operator -D∇ 2 -µ(x) in R N with L-periodicity conditions, where µ(x) = µ + in C + , µ(x) = µ -in C -= C\C + , and µ is extended in R N by L-periodicity.

Integrating equation (2.5) over the cell C and using the periodicity of Φ, we get

- C µ(x)Φ(x)dx = λ 1 [C + , µ + , µ -] C Φ(x)dx. (4.27) Since µ ≥ µ -, it follows that λ 1 [C + , µ + , µ -] ≤ -µ -.
Next, multiply (2.5) by Φ and integrate by parts over C, we get

C D|∇Φ| 2 (x)dx - C µ(x)Φ 2 (x)dx = λ 1 [C + , µ + , µ -] C Φ 2 (x)dx,
and since µ ≤ µ + , we get λ 1 [C + , µ + , µ -] ≥ -µ + . Finally,

-µ + ≤ λ 1 [C + , µ + , µ -] ≤ -µ -.
Step 1: limit of λ 1 [C + , 0, -n] when C + is fixed and n → +∞. Let C + be a fixed Borel subset of C • whose N -dimensional Lebesgue measure belongs to the interval (0, |C|). For n ∈ N\{0}, call µ n the L-periodic function such that

µ n (x) = 0 if x ∈ C + , and µ n (x) = -n if x ∈ C -= C\C + .
Call λ 1,D [C + ] the quantity defined by for all n ∈ N\{0}. Call Φ n the principal eigenfunction of the operator -D∇ 2µ n with L-periodicity, with the normalisation Φ n L 2 (C) = 1. One has

λ 1,D [C + ] = min ϕ∈H 1 0,per (C + )\{0} C + D|∇ϕ| 2 C + ϕ 2 , ( 4
λ 1 [C + , 0, -n] = C D|∇Φ n | 2 + n C - Φ 2 n . (4.30)
From (4.29) and (4.30), it follows that the sequence ( ∇Φ n L 2 (C) ) n∈N\{0} is bounded and that Φ n L 2 (C -) → 0 as n → +∞. Therefore, there exist a subsequence (Φ n ′ ) and a function ϕ ∈ H 1 0,per (C + ) such that Φ n ′ → Φ as n ′ → +∞ strongly in L 2 loc (R N ) and weakly in

H 1 loc (R N ). In particular, Φ L 2 (C + ) = Φ L 2 (C) = 1. Furthermore, since λ 1 [C + , 0, -n ′ ] ≥ D ∇Φ n ′ 2 L 2 (C) , it follows that lim inf n ′ →+∞ λ 1 [C + , 0, -n ′ ] ≥ D ∇Φ 2 L 2 (C) = D ∇Φ 2 L 2 (C + ) .
As a consequence, lim inf

n ′ →+∞ λ 1 [C + , 0, -n ′ ] ≥ λ 1,D [C + ] from (4.28), whence λ 1 [C + , 0, -n ′ ] → λ 1,D [C + ]
as n ′ → +∞, from (4.29). Eventually, since the limit is uniquely determined, one concludes that the whole sequence converges, that is

λ 1 [C + , 0, -n] → λ 1,D [C + ] as n → +∞.
Step 2: A pointwise inequality for the principal eigenfunction associated to an optimal configuration. Let µ -< µ + and h ∈ (0, 1) be fixed. From Appendix A, there exists a Borel set C + ⊂ C • such that C + ∈ S h,µ + ,µ -(S h,µ + ,µ -is defined by (2.9)). Notice that the set C + could be chosen to be symmetric with respect to the hyperplanes {x i = L i /2} for all 1 ≤ i ≤ N , but this property is not needed here. Call µ the L-periodic function which is equal to µ + in C + and to µ -in C -= C\C + . Let Φ be the (unique up to multiplication) principal eigenfunction of the operator -D∇ 2 -µ with L-periodicity. The function Φ is positive, it belongs to all Sobolev spaces W 2,p loc (R N ) for 1 ≤ p < +∞ (in particular, it is of class C 1,α (R N ) for all 0 ≤ α < 1) and it satisfies

-D∇ 2 Φ -µ(x)Φ = λ 1 [C + , µ + , µ -]Φ.
The goal of this step is to prove that inf

(C + ) • Φ = sup (C -) • Φ. (4.31) It is enough to prove that Φ(x 0 ) ≥ Φ(y 0 ) (4.32)
for all x 0 and y 0 which are interior points of C + and C -respectively. There exists ε 0 > 0 such that B ε (x 0 ) (resp. B ε (y 0 )) is included in the interior of C + (resp. in the interior of C -) for all 0 < ε < ε 0 , where B r (x) denotes the open Euclidean ball of centre x and radius r > 0. For ε ∈ (0, ε 0 ), call

C + ε = C + \B ε (x 0 ) ∪ B ε (y 0 ) and C - ε = C\C + ε .
By construction, the Borel set C + ε is included in C • and its Lebesgue measure is equal to h|C|. Let µ ε be the L-periodic function which is equal to µ + in C + ε and to µ -in C - ε , and let Φ ε be the principal eigenfunction of -D∇ 2 -µ ε with L-

periodicity. Call λ = λ 1 [C + , µ + , µ -] and λ ε = λ 1 [C + ε , µ + , µ -] (for 0 < ε < ε 0 ). By definition of λ, there holds λ ≤ λ ε . (4.33)
The functions Φ and Φ ε are positive and L-periodic solutions of

-D∇ 2 Φ -µ(x)Φ = λΦ -D∇ 2 Φ ε -µ ε (x)Φ ε = λ ε Φ ε in R N (4.34)
and up to normalisation, one can assume that Φ ∞ = Φ ε ∞ = 1. The families (λ ε ) 0<ε<ε0 and (µ ε ) 0<ε<ε0 are bounded in R and L ∞ (R N ) respectively. From standard elliptic estimates, the family (Φ ε ) 0<ε<ε0 is then bounded in all W 2,p loc (R N ) spaces with 1 ≤ p < +∞. Standard arguments (see for instance [START_REF] Berestycki | Analysis of the periodically fragmented environment model : I -Species persistence[END_REF]) imply then that λ ε → λ and Φ ε → Φ in C 1 (R N ) (at least) as ε → 0. Multiply now the first equation of (4.34) by Φ ε and the second one by Φ and subtract and integrate over C the two quantities. It follows that

- C (µ -µ ε )ΦΦ ε = (λ -λ ε ) C ΦΦ ε ≤ 0 from (4.

33) and the positivity of Φ and Φ

ε . In C, the function µ -µ ε is equal to µ + -µ -in B ε (x 0 ), to µ --µ + in B ε (y 0 ) and to 0 otherwise. Since µ + > µ -, it follows that Bε(x0) ΦΦ ε ≥ Bε(y0) ΦΦ ε .
By dividing the previous inequality by α N ε N (where α N denotes the Lebesgue measure of the unit Euclidean ball in R N ) and passing to the limit as ε → 0, one concludes that Φ(x 0 ) 2 ≥ Φ(y 0 ) 2 . Since Φ is positive, the inequality (4.32) follows and the claim (4.31) is proved.

Notice that the inequality (4.31) implies in particular that the function Φ is constant on ∂(C + ) • ∩ ∂(C -) • . This property was observed in the numerical computations which were carried out in Section 3.

Step 3: Conclusion. Up to rotation of the frame, one can assume that L 1 = min(L 1 , . . . , L N ). Call r = L 1 /2 and c = (L 1 /2, . . . , L N /2) the centre of the cell C. The ball B r (c) is included in C and it is tangent to the two hyperplanes {x 1 = 0} and {x 1 = L 1 }, like the black region in Fig. 6. Call α N r N the Lebesgue measure of this ball, and assume that

h = α N r N |C| .
Assume now that there is a Borel set C + ⊂ C • such that C + ∈ S h,0,-n for all n ∈ N\{0}. In other words, we assume that

λ 1 [C + , 0, -n] = min K∈S h λ 1 [K, 0, -n] for all n ∈ N\{0}. (4.35)
From this assumption, we claim that there exist a point x ∈ C and a zeromeasure set E such that

C + ⊂ B r (x) ∪ E.
(Φ * ) 2 < C + D|∇Φ| 2 C + Φ 2 = λ 1,D [C + ]
since C + is not included in any ball of radius r, up to a zero-measure set. Notice in particular that the restriction of Φ

* to B r (c) is in H 1 0 (B r (c)). Now define Φ * = Φ * in B r (c) and Φ * = 0 in C\B r (c), and extend Φ * in R N by L-periodicity. The function Φ * is then in H 1 0,per (B r (c))\{0}, whence λ 1,D [B r (c)] ≤ Br(c) D|∇Φ * | 2 Br(c) (Φ * ) 2 = Br(c) D|∇Φ * | 2 Br(c) (Φ * ) 2 < λ 1,D [C + ],
which yields (4.36).

Step 1 yields λ 1 [C + , 0, -n] > λ 1 [B r (c), 0, -n] for n large enough, which is in contradiction with (4.35), since B r (c) is included in C and has the same measure h|C| as C + . As a consequence, the set C + is included in a ball of radius r, up to a set of zero measure. Up to a shift and a modification with a zero-measure set, one can then assume without loss of generality that where a = (0, L 2 /2, . . . , L N /2) ∈ ∂B c (r). Therefore, the partial derivative of Φ with respect to x 1 at the point a is positive. But since the function µ is symmetric with respect to x 1 in R N (because r = L 1 /2, see Fig. 6), the function Φ is also symmetric with respect to x 1 , by uniqueness. Therefore, ∂Φ ∂x1 (a) = 0, a contradiction has been reached.

C + = B r (c).
As a conclusion, the assumption (4.35) cannot hold, which means that the optimal habitat configurations depend on the difference between the growth rates in habitat and non-habitat regions in general.

Notice that, under the notations in Step 3, the above arguments imply that the ball B r (c) is never in S h,µ + ,µ -, for all µ + > µ -, when h = αN r N |C| . In particular, this shows rigourously that ball-shaped configurations are not always optimal for species persistence.

Remark 4.1 In the bounded case, under the assumptions of Proposition 3.2, if we assume that λ 1 [C + ] = λ 1 , it follows from Step 2 that Φ(x) > Φ(x 0 ) for all x ∈ B • . Then, as in Step 3, applying Hopf's lemma at the point x 0 , we get that ∂Φ ∂n (x 0 ) < 0. This is in contradiction with the boundary conditions.

F. Stripe-shaped configurations are not always optimal

We prove here that stripe-shaped configurations are not always optimal for species persistence. The proof is done in the periodic case, in dimension N = 2, but it can be adapted immediately in higher dimensions. Fix two real numbers µ -< µ + . Call λ D (B R (x)) the principal eigenvalue of -D∇ 2 in the ball B R (x) with Dirichlet boundary conditions on ∂B R (x). Since λ D (B R (x)) = R -2 λ D (B 1 (0)), we can fix R > 0 large enough such that λ D (B R (x)) < µ + -µ - 2 , for all x ∈ R 2 . (4.37)

Call A + = πR 2 , fix two positive real numbers l 1 and l 2 , and a positive real number β 0 such that β 2 0 l 1 l 2 > A + . We will compare the two stripe-configurations which are depicted in Fig. 3, where L 1 = βl 1 and L 2 = βl 2 , with the one for which the habitat is a disc a radius R in the cell C = C β = [0, L 1 ] × [0, L 2 ], for β large. The two stripe-configurations and the disc-shaped configurations are constructed so that the habitat be of fixed area A + in C. We will prove that the disc-shaped habitat is better when β is large.

Pick any β ≥ β 0 , and let x be a point such that B R (x) ⊂ [0, L 1 ] × [0, L 2 ]. From [START_REF] Berestycki | The principal eigenvalue and maximum principle for second order elliptic operators in general domain[END_REF], it is known that λ 1 [B R (x), µ + , µ -] is less than the principal eigenvalue of the operator -D∇2 -µ in the ball B R (x) with Dirichlet boundary condition on ∂B R (x). Since -µ = -µ + in the habitat B R (x), it follows that where the function µ β is L 1 -periodic and, up to shift, one can assume without loss of generality that

µ β (x 1 ) = µ + if 0 ≤ x 1 ≤ A + /L 2 = A + /(βl 2 ), µ β (x 1 ) = µ -if A + /(βl 2 ) < x 1 < L 1 = βl 1 .
The function ψ β is in all spaces W 2,p loc (R) for 1 ≤ p < +∞ and one can assume, up to normalisation, that ψ β ∞ = 1. Furthermore, ψ β is analytic in the open intervals (0, A + /(βl 2 )) and (A + /(βl 2 ), L 1 ). If ψ β reaches its maximum at a point x 0 ∈ (A + /(βl 2 ), L 1 ), then -ψ ′′ β (x 0 ) ≥ 0, whence λ β ≥ -µ β (x 0 ) = -µ -, which is impossible since -µ + < λ β < -µ -. Therefore, the maximum of ψ β is reached in the interval [0, A + /(βl 2 )]. 1 Since the families (λ β ) β≥β0 and (µ β ) β≥β0 are bounded in R and L ∞ (R) respectively, the functions (ψ β ) β≥β0 are bounded in all W 2,p loc (R) for 1 ≤ p < +∞. There exist then a sequence β n → +∞ and a function ψ ∞ which is in all W 2,p loc (R) for 1 ≤ p < +∞, such that λ βn → λ ∞ ∈ [-µ + , -µ -] and the functions ψ βn converge (at least) in C 1 loc (R) to ψ ∞ as n → +∞. Furthermore, since µ βn (x) → µ -almost everywhere in R (the convergence holds for all x = 0), the function ψ ∞ satisfies

-Dψ ′′ ∞ -µ -ψ ∞ = λ ∞ ψ ∞ almost everywhere in R.
The function ψ ∞ is then of class C 2 (R). Furthermore, ψ ∞ (0) = 1 = max R ψ ∞ , whence λ ∞ ≥ -µ -. As a consequence, λ ∞ = -µ -and by uniqueness of the limit, the whole family (λ β = λ V [βl 1 , βl 2 ]) β≥β0 converges to -µ -as β → +∞. Similarly, λ H [βl 1 , βl 2 ] → -µ -as β → +∞. Together with (4.38), we conclude that

λ 1 [B R (x), µ + , µ -] < min (λ V [βl 1 , βl 2 ], λ H [βl 1 , βl 2 ])
when β is large. This means that the two stripe-configurations are not optimal for species persistence when β is large and A + is fixed.

  [START_REF] Pearson | Mathematical Contributions of the Theory of Evolution, XV, A Mathematical Theory of Random Migration[END_REF], n = 15 × 15 = 225, and n + varies between 0 and 225, to investigate all the levels of habitat proportion. The number N conf varies between 1 (n + = 0, 1, 224, 225) to 2, 527, 074 (for n + = 112, 113; compare with C 112 225 ≃ 2.8 • 10 66 ).

Figure 1 :

 1 Figure 1: An example of a configuration C + ∈ S h , with h = 0.2, in the particular case of a domain C = [0, 1] 2 , covered by n = 225 disjoint rectangles, and where C + is the union of n + = 45 of these rectangles. The habitats regions are depicted in black.

p

  be the set obtained by reflecting C +, * b with respect to each axis:

Figure 2 :Figure 3 :

 23 Figure 2: The optimal habitat configuration in terms of the number n + of habitat cells in the case of a bounded environment with Neumann boundary conditions. The black cells correspond to habitat and the white cells to nonhabitat.

Figure 7 :

 7 Figure 7: Comparison between the quarter disc configuration and the striped configuration, in the bounded case of a square C = [0, 1] × [0, 1], in terms of the amplitude of the difference between the growth rates b, and of the habitat abundance h = |C + |/|C|. In the shadowed area, the quarter disc configuration provides lower values of λ 1 , while in white zone the striped configuration gives lower λ 1 values.

  µ * is extended by L-periodicity in R N . From the construction of the functions µ * n , the function µ * satisfies the same symmetry properties as µ * n , in the sense that if x = (x 1 , . . . , x N ) ∈ Ω + , then µ * (y) = µ + for all y ∈ [-|x 1 |, |x 1 |] × • • • × [-|x N |, |x N |]. Observe also that the boundary of Ω + has a zero N -dimensional Lebesgue measure. Let x = (x 1 , . . . , x N ) be an interior point of C \ Ω + . There exists y = (y 1 , . . . , y N ) another interior point in this set, such that y ∈ Q N and |y i | < |x i | for all i = 1, . . . , N. (4.15) Since x, y ∈ Ω + , µ * (x) = µ * (y) = µ -. Moreover, by pointwise convergence, µ * n (y) = µ -for n large enough. From the symmetry properties of µ * n , and from (4.15), we have µ * n (y) ≥ µ * n (x), therefore, µ * n (x) = µ -(= µ * (x)) for n large enough. If now x = (x 1 , . . . , x N ) is an interior point of Ω + , there exists another point y in this set such that y ∈ Q N and |y i | > |x i | for all i = 1, . . . , N , and we similarly obtain µ * n (x) = µ + (= µ * (x)) for n large enough. The above observations imply that µ * n (x) → µ * (x) as n → +∞ almost everywhere in C. Since the functions µ * n are uniformly bounded, Lebesgue's dominated convergence theorem yields µ * n → µ * in L 1 (C) as n → +∞. Therefore, µ * ∈ P h .

  and fix the habitat proportion h = n + n . Let P h be defined by (4.14), and set C + (σ) := {x ∈ C, σ(x) = µ + }, for σ ∈ P h .

  P h be defined by(4.14), and let Ph be the set of functions defined by (4.14) but with C instead of C. For a function σ in Ph , we denote by C+ (σ) := {x ∈ C, σ(x) = µ + }. Let λ 1,b (µ) denote the Neumann principal eigenvalue of -D∇ 2 -µ in C, and λ 1,p (σ) denote the periodic principal eigenvalue of -D∇ 2 -σ in C.

  p in Ph . Indeed, let us set λ 1,b := min σ∈P h λ 1,b (σ) and λ 1,p := min σ∈ Ph λ 1,p (σ), and

  [START_REF] Cantrell | Diffusive logistic equations with indefinite weights: population models in disrupted environments[END_REF] and ψ is L 1 -periodic. Now, set μ(y) := μ(yL 1 ). The function μ is 1-periodic and, in the period cell [0, 1], μ is equal to µ + on an interval of length l′ 1 = l 1 /L 1 = |C + |/(L 1 L 2 ) = |C + |/|C| ∈ (0,[START_REF]IUCN Red List of Threatened Species. A Global Species Assessment[END_REF], and it is equal to µ -on an interval of length 1 -l ′ 1 . The function μ does not depend on L 1 . The 1-periodic function ξ defined by ξ(y) = ψ(yL 1 ) satisfies -(L 1 ) -2 Dξ ′′ (y) -μ(y)ξ(y) = λ V ξ(y), y ∈ R, ξ > 0 in R.(4.24)

  (C + ) denotes the set of H 1 per functions which vanish almost everywhere in C -. When C + is a smooth open set, λ 1,D is the principal eigenvalue of the operator -D∇ 2 in C + with Dirichlet boundary condition on ∂C + ∩ C and L-periodicity. Since µ n = 0 in C + , it follows from (4.28) and (4.26) with µ = µ n that λ 1 [C + , 0, -n] ≤ λ 1,D [C + ] (4.29)

  If not, since the ball B r (c) has the same measure as C + and is included in C, it follows from the definition of λ 1,D [C + ] in Step 1 and from Schwarz rearrangement inequalities that λ 1,D [C + ] > λ 1,D [B r (c)]. (4.36) Indeed, let Φ be a minimiser of (4.28), call Φ its restriction in C and call B R (c) the open ball of centre c and radius R > r which has the same measure as C. The Schwarz decreasing symmetric rearrangement Φ * of Φ with respect to the point c is in H 1 0 (B R (c)) and it is such that Φ * = 0 in B R (c)\B r (c) and Br(c) D|∇Φ * | 2 Br(c)

Fix now any n 0 ∈

 0 N\{0}. Call µ the L-periodic function which is equal to 0 in B r (c) and to -n 0 in C\B r (c), and call Φ the principal eigenfunction of the operator -D∇ 2 -µ with L-periodicity.Call λ = λ 1 [B r (c), 0, -n 0 ] the principal eigenvalue. The function Φ is (at least) of class C 1 (R N ) and it is positive in R N . FromStep 2, the function Φ is equal to a constant β > 0 on ∂B r (c), and it is larger than or equal to β in B r (c). In B r (c), the function Φ solves -D∇ 2 Φ = λΦ and the functionψ = Φ -β satisfies    -D∇ 2 ψ -λψ = λβ > 0 in B r (c), ψ ≥ 0 in B r (c), ψ = 0 on ∂B r (c)since 0 < λ (< n 0 ). It follows from the strong maximum principle and Hopf lemma that ψ > 0 in B r (c) and ∂ψ ∂x 1 (a) > 0.

λ 1 [

 1 B R (x), µ + , µ -] < λ D (B R (x)) -µ + < -µ + + µ - 2 (4.38)from our choice of R in (4.37).Let us now estimateλ β := λ V [L 1 , L 2 ] = λ V [βl 1 , βl 2 ], for β ≥ β 0 .From the arguments used in Appendix E, there exists a positive function ψ β defined in R, which is L 1 -periodic and satisfies-Dψ ′′ β -µ β (x 1 )ψ β = λ β ψ β almost everywhere in R

Actually, this property could be viewed as a consequence of Step

of Appendix D, since µ β corresponds to an optimal configuration among all L 1 -periodic one-dimensional configurations for which the habitat has a total length equal to A + /(βl 2 ) in [0, L 1 ].

Acknowledgements

The authors are grateful to Etienne K. Klein for his valuable comments and suggestions.

Appendices

This section is devoted to the proofs of the mathematical results which were stated throughout this paper. The results of the appendices A,B,D,E,F are true in both periodic case and bounded case when Ω is a rectangle. However, we decided, for the sake of simplicity, and since the proof in one case can easily be adapted to the other, to prove each result in only one case.

A. Existence of an optimal configuration

The existence of at least one optimal configuration that maximises the chances of species persistence is proved in the periodic case. For another proof, in the Dirichlet case, see [START_REF] Chanillo | Symmetry Breaking and Other Phenomena in the Optimization of Eigenvalues for Composite Membranes[END_REF].

Let µ -< µ + be two given real numbers, L 1 , . . . , L N be N given positive real numbers, and h ∈ [0, 1]. Call L = (L 1 , . . . , L N ) and


Gilles Dodinet
email: gilles.dodinet@lamsade.dauphine.fr

Michel Zam
email: zam@lamsade.dauphine.fr

Genevieve Jomier
email: genevieve.jomier@lamsade.dauphine.fr

Geneviève Jomier

COEVOLUTIVE META-EXECUTION SUPPORT Towards a design and execution continuum

Keywords: Software evolution, collaborative development, model-driven engineering

  L'archive ouverte pluridisciplinaire

INTRODUCTION

Model Driven Engineering (MDE) considers models as first-class citizens in the development process. However, consistent evolution, and thus coevolution, of target systems, their models and their meta-models is both a necessity and a challenge.

The last decade has seen a rapid expansion of MDE techniques and applications, partly due to hopes of knowledge capitalization. While software development process is a "multilevel, multiloop feedback system" involving a variety of actors [START_REF] Lehman | Software's Future: Managing Evolution[END_REF], the main responsibility still relays on the developer's shoulders. Although the stakeholder has the best understanding of the functional requirements, he is still kept away from the whole development process.

In order to reduce the effect of software aging [START_REF] Grubb | Software Maintenance Concepts and Practices[END_REF], and the associated maintenance cost (about 50% to 75% of the global development cost, McKee, 1980), the stakeholder must be actively involved in the software process [START_REF] Ambler | Agile Modeling: Effective Practices for Extreme Programming and the Unified Process[END_REF]. As a consequence, in a MDE context, any model evolution should be immediately reflected in the system to assert its validity. To support such requirements, ad hoc tools must be provided to merge the development and execution phases. [START_REF] References Agrawal | The Claremont report on database research[END_REF] acknowledges the need for tools that would give greater power to stakeholder while empowering the developer.

The coevolution of models with their own models (called meta-models) and with their instances remains a critical issue. (Van Deursen et al., 2007) points the following paradox: while MDE approaches have been designed to handle continuous evolution there's today little support for model and meta-model coevolution. This paradox also applies to coevolution of models and their instances -i.e. the systems represented by these models. Today models are essentially static artifacts: they usually do not evolve easily after the initial development phase. Once a system is deployed, any modification made to the structure of its model cannot be propagated to the running system without going through another full development cycle.

This article presents a shared execution environment built upon a collaborative modeling support that intends to answer the general problem of coevolution of models and their instances. The system is decomposed into fine-grain elements common to the various abstraction levels. They are used to express both models and their instances, including application data and programs.

The article focuses on data-centered applications using relational database to represent models and meta-models stored as user data and metadata.

The next section describes a typical evolution scenario spanning all abstraction layers; section 3 presents related work; section 4 presents the proposed solution and section 5 describes an implementation. Section 6 concludes and opens perspectives.

COEVOLUTION SCENARIO

As requirements change, modifications can occur at any abstraction level and must be propagated consistently. The simple, yet representative, scenario below illustrates such coevolution. Figure 1 shows the initial state of the scenario introducing the Employee class and an instance, e. Employee is also an instance of the Class class and defines the hello instance method. A call to hello prints the class structure (via a call to the show method defined in the Class class), followed by the attribute values of e. To simplify, the Attribute class does not appear. UML 2.x class diagram notation is used and the pseudo-code is inspired by Python syntax.

M0, M1 and M2 represent abstraction layers as defined by the OMG modeling stack [START_REF] Omg | Architecture Board ORMSC. Model driven architecture (MDA). OMG document number ormsc[END_REF]. M2 is the meta-model layer. M1 is the model layer. Any model element defined in M1 is an instance of a meta-element defined in M2. M0 is the system layer; every system element in M0 is an instance of a model element defined in M1.

Layers M3 (meta-meta-model level) and M2 are merged.

The goal is to be able to modify all three layers dynamically, while the system is still running, and take the modifications into account immediately without restarting or redeploying the system.

Figure 2 shows the result of important evolutions: the introduction of the single inheritance concept used by the Employee class, while its instance is subject to an independent update (salary raise). This is achieved by executing the following operations: The result of the execution of e.hello() is shown on the right hand of the M0 layer, Figure 2.

This simple example demonstrates that modifications can affect all the abstraction layers of a system. It shows how flexible a coevolutive system needs to be. More complex scenarios may require a developer decision regarding whether and how to propagate the changes through the abstraction layers.

Such a scenario states the requirement of model executability. It underlines the necessity of coexistence of target data with models and metamodels, and the need for merging design-time and runtime.

RELATED WORK

This section presents existing work regarding three main areas related to the coevolutive support.

Model executability

Although it is not a recent research trend, model and meta-model executability is getting more and more attention [START_REF] Mellor | Executable Uml: a Foundation for Model-Driven Architectures[END_REF]. (Breton and Bezivin, 2001) attempts to identify the main characteristics of executable meta-models, and underlines the need of an additional layer to specify the operational semantics of a model.

Magritte [START_REF] Renggli | Magritte -A Meta-Driven Approach to Empower Developers and End Users[END_REF]) is a reflexive meta-description framework refining SmallTalk reflexive It allows describing the business classes of an application, as well as their respective attributes, operations and constraints. It abstracts the Type-Object design pattern [START_REF] Johnson | Type object[END_REF] and blurs the boundary between component and property concepts. Magritte primarily targets web applications and dynamically interprets the metadata to generate the user interface. [START_REF] Ducasse | Meta-Environment and Executable Meta-Language using Smalltalk: an Experience Report[END_REF]) considers metamodel executability and expands the applicability of Magritte's meta-description to embrace system classes.

While promising, these approaches based on source code, are developer oriented. Other approaches include the use of UML profiles and static model transformations [START_REF] Hemel | Code Generation by Model Transformation[END_REF] to generate fully functional applications or prototypes. However, the approaches presented here do not cover data and model consistent coevolution.

Schema and data management

Data consistency is a property guaranteed by any DBMS, stating that database changes, called transactions, are processed reliably. However, traditional database schema evolution remains complex and heavy and do not match with the fine grained and frequent rhythm of MDE modeling evolutions. [START_REF] Dinu | Guidelines for the effective use of entity-attribute-value modeling for biomedical databases[END_REF] Using this approach, a model modification does not require the update of the underlying database schema. The model can then easily be adapted (as far as only data are concerned). EAV is often used when data structure must change frequently or is not known in advance.

Although the EAV approach allows modifying easily the virtual target schema, the absence of even basic model constraints is a major drawback. Therefore, EAV with Class and Relationships model (EAV/CR, [START_REF] Nadkarni | Organization of heterogeneous scientific data using the EAV/CR representation[END_REF] has been introduced. It refines the EAV approach with schema concepts like class and relationship. Specific tables are used to store metadata. This allows refining the schema dynamically. However this solution is not reflexive. As a consequence, metadata cannot evolve without updating the relational structure.

It is also important to note that those approaches (EAV and EAV/CR) do not include operational semantics: they only deal with persistent data structure, while the behavior of stored objects is not defined.

Software evolution management

Source Configuration Management (SCM) tools, such as CVS, are mainly used to trace software evolutions. They are well suited for general textbased artifacts (e.g. java or xml files). For instance they can be used, in conjunction with bug trackers, to trace evolution and get a better understanding of the past changes [START_REF] D'ambros | Analyzing software repositories to understand software evolution[END_REF].

From a MDE perspective, coevolution refers to the need for the models to evolve with their metamodels. However, as files size is usually much bigger than the size of model elements, SCM are not well adapted to the fine grain management of the evolution of strong typed model elements and change impact. (Wacshmuth, 2007) proposes a classification of meta-models based on semantics preservation properties. Models are incrementally adapted through a serie of changes represented by highlevel transformations. Those transformations are executed manually, making the successive evolution steps explicit. [START_REF] Hôßler | Coevolution of models, metamodels and transformations[END_REF] presents a set of usual transformations used at various abstraction levels to adapt a meta-model and to migrate its models automatically. Meta-models are supposed to be instances of a reference meta-meta-model, i.e. the M3 layer in the OMG modeling stack. Every model element is connected to its previous version using a predecessor relationship.

None of these solutions consider model executability nor target application data.

The Karma model [START_REF] Zamfiroiu | La traçabilité du processus de conception en génie logiciel[END_REF] integrates typical software configuration features (such as check-in, tags, branches, etc.) into relational databases. Karma is based on a version model independent of the data model, called DBV, that was first introduced in [START_REF] Cellary | Consistency of Versions in Object-Oriented Databases[END_REF]. Any data modification is automatically traced and any previous state of the system can be restored consistently.

The solutions presented in this section only address one particular facet of the actual problem. In order to handle evolution at any abstraction level throughout the whole software lifecycle, all three major aspects should be addressed conjointly in a consistent and generic way: model executability, schema and data management, and finally consistent software evolution management.

COEVOLUTION SUPPORT

This section introduces a design and execution platform integrating a coevolution support. All abstraction layers coexist within the runtime environment. Every system element is decomposed into elementary constructs, embracing the EAV strategy. The coevolution of those elements is enabled by the Karma model.

We successively cover the coevolution underlying model, the persistency management and the traceability of evolutions, designed to fulfill the evolution requirements introduced in section 2.

Coevolution model

In order to understand how to build the integrated coevolution model, a semantics-free descriptive layer (atomic level) is first introduced. Above it stays a logically typed structural semantics level, called molecular. Finally, an operational layer is added dynamically.

Atomic and molecular layers

The descriptive layer sets the foundations of the coevolutive execution model. It provides the building blocks of the system. By convention, any element of the system is described using undividable units, called Atoms. Taken separately atoms do not have additional meaning. The business semantics will only emerge from their values and their composition. We call this model "atomic model" by analogy with the ancient Greek atomic model (Figure 3). An Atom represents (at most) one element of the system and holds its identity. It always references a meta Atom, which holds the identity of the element describing the structure and behavior of its associated element, even if at this stage we don't know how to interpret the structure. An Atom has a collection of Slots. A Slot represents a primitive valued property (string, integer, etc…). The Spin is a particular type of Slot. Its value is a reference to an Atom.

Pursuing the atomic metaphor, Atoms can be composed arbitrarily to build more complex constructs, called Molecules. An Instance is a particular molecule that represents an object in the Class-Object paradigm.

In order to be correctly interpreted and validated (structurally and behaviorally), molecules must be typed.

Let us consider the Figure 4 below, representing a simple view on the M2 abstraction layer -where Class and Attribute are classes. Could they be represented using atoms and molecules?

Operational semantics layer

In order to make the model executable the Class and Attribute classes should be equipped with operational semantics. The operational semantics should express how a particular atomic assemblage must be interpreted and transformed into a programmatic molecular construct, e.g. how to reconstruct the Figure 4 from Figure 5.

Thus, the meta relation is equipped with an instantiation semantics defined by the simple rules below (Table 1). T is the transformation operation that creates a molecule from an atom. A is the set of atoms present in the system. instanceof is an operator defined at the molecular level, returning true if the second operand is a representation of the first operand in the upper abstraction level. The operational semantics must then assert that the above rules are always respected. These constraints can be enforced when an instance of the Molecule class is created (see Figure 3), or later through explicit validation.

(1) ∀ a2 ∈ A ∃! a1 ∈ A : a2.meta=a1 (2) ∀ a1, a2 ∈ A
Although an additional layer is needed to specify this operational semantics, it can be added afterwards. Thus, meta-modeling activity can be decomposed into two distinct phases: (i) assemble the atoms, giving them a form that reflects the foreseen semantics, (ii) give the assemblages an explicit semantics and behavior through the operational layer that can be refined later dynamically.

If atomic structural requirements change, the operational semantics must be revised as well to consider the structural modifications. This is possible because the definition of the operational layer is dynamic. Both operational semantics specification and specific behavior implementation are represented as model operations and as such decomposed into Atoms as well.

Also, since all abstraction layers have an atomic representation, they all coexist. So the system as a whole becomes causally connected, i.e. any change made to the self-representation of the system is immediately reflected in its actual state and behaviour, as defined in [START_REF] Maes | Computational reflection[END_REF].

Persistency and evolution

We are finally able to consider the persistency of atoms and molecules, and how their coevolution is supported.

As every element of the system is ultimately represented as an atom, it is possible to use a uniform solution to store the meta-meta-model, all meta-models, all models, and all terminal instances as mere data in a relational database.

For instance, atoms, as well as slots and spins (as described in Figure 3) can be stored in a unique table. In any case, this database has a fixed schema. This remarkable property eliminates common problems due to schema evolution and data migration. As a consequence, restarting and redeploying the system is no longer mandatory.

Moreover, querying and manipulating data and models become uniform and can be expressed at the atomic level. Thus, specific semantics must be used in order to interpret the atoms. For instance, since there is no Employee table anymore, the evaluation of a query like select * from Employee must be adapted. Therefore, the system can run continuously while still evolving.

Tracing the evolution

The database is extended with a traceability mechanism, based on the Karma and DBV models (Section 3.3). Thus not only the current state but also all the previous states of the database are available in conjunction with the user operation traces. The database states -successive as well as alternative -are organized as a global version tree (Cellary and Jomier, 2000). The database is therefore multi-version. Each version of the database, called global version, can be considered as a consistent mono-version database and contains at most one version of each atom. In the Karma traceability model, any elementary evolution creates a new system state automatically, derived from the previous state. Each local version of a given Atom is marked with the global version identifier where the evolution occurs. Therefore, it becomes easy to extract any previous state using a version identifier.

It is also possible to derive explicitly an alternative branch from a given global version. Finally, any given previous state of the database can be recovered instantaneously by selecting the corresponding version identifier in the version tree.

Molecules are always assembled from appropriate Atom versions in a given global version of the database. A system state is based on a given database version. Thus, its atoms and molecules remain consistent. Since models are composed of Molecules as well and considered as mere data (and are treated as such), they are also intrinsically versioned. That means that all the abstraction layers -M0, M1, and M2 -benefit in the same way from the multi-version properties.

Some radical evolutions like splitting or merging elements can alter the molecular identity. Even in this case, traces and versions help detecting and avoiding inconsistencies.

Finally, not only the system can run continuously while evolving, but its entire history is recorded and the evolution can follow different branches. A continuum of design and execution lifecycle phases is thus achieved as long as every resource is in fine equipped with an atomic representation. Then the version substitution principle can be applied and the system can improve continuously.

IMPLEMENTATION

General architecture

The atomic model presented above has been implemented as an integrated modeling platform, based on JavaEE. The user interface paradigm follows the single-page web application style, as defined in [START_REF] Mesbah | Migrating Multi-page Web Applications to Single-page Ajax Interfaces[END_REF]. The platform allows the stakeholder and the developer to create models that can be dynamically instantiated. Their meta-models can also be refined dynamically. Any modification is propagated immediately all the way to target data.

The meta-model executability is achieved through a bootstrap phase during which the system operational semantics is dynamically injected, using a script. Two different implementations of the bootstrap script have been written (in Javascript and in Python) and can be run concurrently on the server side (JSR-223). Figure 6 sketches the environment architecture and interactions. Let us suppose the user wants to execute a simple script that instantiates the Employee class. Users interact with the application using standard internet browsers and specify the operation to be executed remotely (as a script). First, the bootstrap script is evaluated to create the required execution context. Then, the user action is executed on the server side. This is necessary in order to transform a Molecule into a python or javascript object through a SemanticsInterpreter component and to address it in a way that makes functional sense. SemanticsInterpreter delegates the Molecule construction to the InstanceManager component. Therefore, from a developer perspective, the Molecule structures are hidden behind scripted proxies, as illustrated in section 5.2.

Every scripting object that has an underlying molecular description (section 4.1.1), wraps its corresponding Molecule (Java) instance; any read or write operation applied on such a scripting object is redirected to its Java counterpart. This also allows adding strong typing to usually loose typed languages (i.e. Javascript or Python).

The Karma model has been implemented at two different levels.

-In the database: additional tables have been created to store the global versions and the branches; the Atom table has been changed to a view in order to generate Atom versions transparently. Finally a set of stored procedures has been introduced to navigate through the timeline and the global version tree.

-In the application framework: an API has been added at the Java side and in the bootstrap script to enable global version selection.

Application specific behavior, as well as metamodel behavior, is implemented through class operations, which are instances of the Operation class. Thus, operations are also considered as mere data. Their body is stored in a Text column. Their dependencies can then be detected both statically, by parsing their body, and dynamically, by automatically collecting their execution traces, explicitly managed by the Karma model.

Example

The solution implemented on top of the proposed model covers the scenario presented in Section 2, as illustrated in Figure 7 and Figure 8 -those are screenshots of a scripting console embedded in the actual platform. First the Employee class and its instance called popeye are created and made persistent. At this stage, the hello method is called on the instance popeye. The result is shown on Figure 7. Then the Person class is created and the Employee class is redefined as a specialization of Person (by setting its "super" attribute). Next, the nickname attribute is moved under the Person class and the show method is dynamically updated. Finally the instance is loaded again and its salary is raised. Figure 8 shows the result of calling again the hello method on the Popeye instance.

CONCLUSION

This article introduces a coevolution support for reflective meta-models and their instances, tightly integrated into an MDE execution platform. Any system element is decomposed into undividable units, called atoms, allowing a homogeneous representation of any element, regardless its abstraction level.

To illustrate our approach we chose to reproduce the Class-Object paradigm. While the Class-Object approaches mainly stand at the conceptual level, Atoms and Molecules constructs are introduced at a lower abstraction level to allow refining (meta-)models without having to alter the underlying storage structure.

Atomic representation of both models and their instances are stored in a shared multi-version database. Since models and instances coexist as high-level atom constructs, their coevolution is managed in a consistent way, and model evolutions are instantly reflected into their instances, making the whole system causally connected.

The prototype we implemented allows stakeholders, developers and final users to define, update and run models and theirs instances concurrently. It has been experimented in a multi criteria decision aid (MCDA) platform called DECISIONDECK. Practitioners design formal MCDA methods using basic user interface to edit and visualize input data. Although more formal evaluation is needed the preliminary results are encouraging.

However, our solution implies a development paradigm shift, and as such requires appropriate tools. We need now to focus on the development environment in order to support the usual industrial constraints of quality and productivity. These enhancements will be implemented in the next version of the system mainly as molecular constructs and will become a part of its karma.

Figure 1 :

 1 Figure 1: Employee class, version 1

-

  Add a reflexive association, super, to the Class class (M2 layer) -Accordingly, update the show method to list the inherited attributes (M2 layer) -Create the Person class generalizing the Employee class (M1 layer) -Move the nickname attribute from Employee to Person (M1 layer) -Update the value of e.salary (M0 layer)

Figure 2 :

 2 Figure 2: Employee class, version 2

Figure 3 :

 3 Figure 3: Atomic model

Figure 4 :

 4 Figure 4: M2 layer -Class and Attribute classes. An atomic and molecular solution is depicted in Figure 5 below. The Class class is an instance of itself, and the attributes of the Attribute class are instances of Attribute. The « instance of » relationship is concretized in the atomic model by a "meta" relation between the Atom representing a given object and the Atom holding the identity of its class. Finally, the Class and Attribute classes are interpreted as typed molecules aggregating necessary atoms. The aggregation between the Class class and the Attribute class has no

Figure 5 :

 5 Figure 5: Atomic and molecular representations of the Class and Attribute classes

  ∀ a1, a2, a3, a4 ∈ A T(a2) instanceof T(a1) and T(a3) instanceof T(a4) and a1.meta=Class and a4.meta=Attribute and a3.owner=a2 and a3.type='Slot' ⇔ T(a2)[a4.name] = a3.value

Figure 6 :

 6 Figure 6: Implementation: architecture overviewSignificant modifications of the meta-model might need to be reflected in the bootstrap script, so the script must remain modifiable dynamically.Figure6sketches the environment architecture and interactions. Let us suppose the user wants to execute a simple script that instantiates the Employee class. Users interact with the application using standard internet browsers and specify the operation to be executed remotely (as a script). First, the bootstrap script is evaluated to create the required execution context. Then, the user action is executed on the server side. This is necessary in order to transform a Molecule into a python or javascript object through a SemanticsInterpreter component and to address it in a way that makes functional sense. SemanticsInterpreter delegates the Molecule construction to the InstanceManager component. Therefore, from a developer perspective, the Molecule structures are hidden behind scripted proxies, as illustrated in section 5.2.

Figure 7 :

 7 Figure 7: Employee instance creation

Figure 8 :

 8 Figure 8: M0-M2 layers coevolution Meta-model evolution is illustrated Figure 9, ci-dessous through the update of the show method of the Class class. This evolution doesn't change how Atom constructs are interpreted and thus doesn't require reconsidering the bootstrap phase (although the introduction of the inheritance concept does).

Figure 9 :

 9 Figure 9: show operation

  propose a more flexible modeling technique called the Entity-Attribute-Value (EAV) model. EAV allows defining and storing application data in a relational database without knowing the target schema beforehand. Application data are stored in three tables: the Entity table holds the object identities; the Attribute table holds the attribute definitions; finally the Value table holds the actual values for a pair Entity/Attribute.

Table 1 :

 1 Instantiation rules