Matthieu Alfaro

Thomas Giletti

Yong-Jung Kim

Gwenaël Peltier

Hyowon Seo

ON THE MODELLING OF SPATIALLY HETEROGENEOUS NONLOCAL DIFFUSION: DECIDING FACTORS AND PREFERENTIAL POSITION OF INDIVIDUALS

Keywords: nonlocal diffusion, heterogeneity, focusing kernels, deciding factors, shape of steady states AMS Subject Classifications: 92B05 (General biology and biomathematics), 45K05 (Integro partial differential equations), 35B36 (Pattern

We develop general heterogeneous nonlocal diffusion models and investigate their connection to local diffusion models by taking a singular limit of focusing kernels. We reveal the link between the two groups of diffusion equations which include both spatial heterogeneity and anisotropy. In particular, we introduce the notion of deciding factors which single out a nonlocal diffusion model and typically consist of the total jump rate and the average jump length. In this framework, we also discuss the dependence of the profile of the steady state solutions on these deciding factors, thus shedding light on the preferential position of individuals.

Introduction

The purpose of this exploratory paper is to investigate the modelling of heterogeneous nonlocal diffusion and the connection with local diffusion. Starting from the most general nonlocal heterogeneous diffusion, we will see how to obtain a large class of local heterogeneous diffusion equations by taking a focusing kernel limit. This process sheds new light on the underlying dynamics of these models. In particular, the notion of deciding factors, and its implications on the shape of the steady state solutions, may help to find out whether a heterogeneous diffusion model is appropriate or not for a given application.

1.1. Homogeneous diffusion. Under a homogeneous environment, a local diffusion equation is often given by the classical heat equation (1.1)

u t = D∆u,
where the diffusivity D > 0 is constant. This diffusion equation is used in a wide range of applications from physics problems to ecology ones, where the homogeneous diffusion, D∆u, describes a random dispersal and x ∈ R N the spatial position. Such a linear diffusion model can be recovered from a probabilistic individual based model. In the study of evolutionary biology, the linear diffusion also models the mutation process, where x ∈ R N denotes a phenotypic trait (see [START_REF] Alfaro | Explicit solutions for replicator-mutator equations: extinction versus acceleration[END_REF][START_REF] Fisher | The wave of advance of advantageous genes[END_REF][START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF][START_REF] Kolmogorov | Etude de l'équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique[END_REF][START_REF] Shigesada | Biological invasions: theory and practice[END_REF][START_REF] Skellam | Random dispersal in theoretical populations[END_REF] among many other references).

Remember that the diffusion equation (1.1) is for an environment where there is no advantageous direction (isotropic) nor location (homogeneous).

In ecology, evolutionary biology, epidemiology, etc. the model usually comes with a nonlinear reaction term, which takes into account birth, death, competition, etc. Throughout this paper, we neglect the nonlinear population dynamics and focus on the role of diffusion.

When the spatial domain is the whole Euclidean space R N , the nonlocal counterpart of the heat equation (1.1) is

(1.2) u t = R N K(x -y)u(t, y) dy -mu(t, x),
where K : R N → R + is a nonnegative integrable kernel and the coefficient m is the total jump rate,

m = R N K(z) dz > 0.
The way to recover the local equation (1.1) from the nonlocal one (1.2) through a focusing kernel limit is classical and will be made precise in Section 3. The integro-differential equation (1.2) is motivated by the long-range dispersal phenomenon meaning that an individual may "travel" long distance in a short time scale. Nonlocal diffusions are therefore very popular in population models for which long-distance dispersal events are relevant. The propagation of epidemics [START_REF] Medlock | Spreading disease: integro-differential equations old and new[END_REF][START_REF] Murray | Mathematical biology II: spatial models and biomedical applications[END_REF] and dispersal of seeds [START_REF] Clark | Why trees migrate so fast: confronting theory with dispersal biology and the paleorecord[END_REF] are such examples. Note that, in evolutionary biology, the mutation is also typically modeled in terms of integral operators [START_REF] Alfaro | Confining integro-differential equations from evolutionary biology: ground states and long time dynamics[END_REF][START_REF] Bürger | On the maintenance of genetic variation: global analysis of Kimura's continuum-of-alleles model[END_REF][START_REF]Perturbations of positive semigroups and applications to population genetics[END_REF][START_REF] Fleming | Equilibrium distributions of continuous polygenic traits[END_REF][START_REF] Gil | Mathematical properties of a class of integro-differential models from population genetics[END_REF][START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF][START_REF] Lande | The maintenance of genetic variability by mutation in a polygenic character with linked loci[END_REF]. In the following, we mainly refer to "spatial dispersal" but, obviously, our analysis also applies to "mutations".

In the context of a spatial dispersal, the kernel K(z) stands for the jump rate for an individual to move from x to x + z or, equivalently, K(x -y) is the jump rate from y to x. Hence, the constant m is the total jump rate from a given point x ∈ R N . In particular, the last term in (1.2) models the amount of species departing the position x. By integrating (1.2), one may check that the total population of the solution is preserved, i.e.

R N u(t, x)dx is constant.

The central underlying assumption in the nonlocal diffusion model (1.2) is that the jump rate from x to x + z is independent from the departure point x and the arrival one x + z. In most mathematical literature, the constant m is normalized to 1 which can be done without loss of generality and the function K is assumed to be radial (see [START_REF] Alfaro | Propagation phenomena in monostable integrodifferential equations: acceleration or not?[END_REF][START_REF] Andreu-Vaillo | Nonlocal diffusion problems[END_REF][START_REF] Chasseigne | Asymptotic behavior for nonlocal diffusion equations[END_REF][START_REF] Coville | On a non-local equation arising in population dynamics[END_REF][START_REF] Garnier | Accelerating solutions in integro-differential equations[END_REF] and references therein). As a result, the obtained nonlocal equation (1.2) is isotropic and homogeneous in space, as is (1.1).

1.2. Heterogeneous diffusion. However, the motility of biological species (in both local and nonlocal cases) depends on the environment, which is usually spatially heterogeneous. This spatial heterogeneity gives us a challenge to develop diffusion models that take into account the spatial heterogeneity properly in local and nonlocal equations. One of our goals is to address both issues jointly and to draw new connections between those models.

To give the motivation, let us first illustrate the challenge in modelling heterogeneous diffusion starting with a local diffusion case. A naive approach would be to simply replace the diffusivity D in (1.1) by a real-valued function D(x) and obtain "u t = D(x)∆u". However, the resulting evolution equation does not satisfy the mass conservation property (even under the zero-flux boundary condition) and one should instead consider a diffusion model in a divergence form for the motion of biological organisms.

On the other hand, there are infinitely many choices of diffusion equations in a divergence form when the diffusivity D is not constant, which can be written

(1.3) u t = ∇ • D(x) q ∇(D(x) 1-q u) , q ∈ R.
Obviously, when the diffusivity coefficient D is constant, these diffusion laws are all equivalent to the heat equation (1.1). In a heterogeneous case, they lead to different solutions. Furthermore, if D is not smooth enough, the "classical weak solution" is not defined and a different approach, depending on the choice of the diffusion law, is required. Many well-known diffusion equations are written in the form of (1.3). First, one recovers the so-called Fick's diffusion law [START_REF] Fick | On liquid diffusion[END_REF] when q = 1:

u t = ∇ • D(x)∇u)), (Fick's law).
This is a case without an advection phenomenon and constants are steady states. When q = 0, the equation is written as

u t = ∆(D(x)u). (Chapman's law)
Chapman derived this diffusion law using the kinetic theory in [START_REF] Chapman | On the brownian displacements and thermal diffusion of grains suspended in a non-uniform fluid[END_REF] where his main purpose is to explain the thermal diffusion phenomenon. Before Chapman's derivation, Wereide [START_REF] Wereide | La diffusion d'une solution dont la concentration et la température sont variables[END_REF] derived a diffusion law that corresponds to the case with q = 1 2 ,

u t = ∇ • D(x)∇( D(x)u) , (Wereide'
s law) using basic physical intuitive relations. The latter two equations are also respectively called Itô and Stratonovich type since they are satisfied by probability density functions of the corresponding stochastic interpretations of the Brownian motion. Unlike Fick's law, the two contain drift phenomena and hence they may also be called Fokker-Planck equations. The diffusion laws in (1.3) satisfy the mass conservation property, and it is not obvious to see which one is the right diffusion law among all possibilities. The related question of whether heterogeneous dispersal is advantegeous to a species has been adressed from the evolutionary point of view in several works [START_REF] Averill | The role of advection in a two-species competition model: a bifurcation approach[END_REF][START_REF] Cantrell | Movement toward better environments and the evolution of rapid diffusion[END_REF][START_REF] Cantrell | Evolution of dispersal in heterogeneous landscapes[END_REF]. Recently, two of the authors of the present work derived heterogeneous diffusion equations from revertible kinetic systems and compared the properties of diffusion laws using a thought experiment [START_REF] Kim | Model for heterogeneous diffusion[END_REF]. In particular, they obtained Wereide's diffusion law as a special case. In this paper we see that the same questions and issues exist and need to be addressed for the spatially heterogeneous nonlocal diffusion case. In this process we hope to obtain new insights on the long standing issue of heterogeneous diffusion (see, e.g., [START_REF] Novak | Habitat heterogeneities versus spatial type frequency variances as driving forces of dispersal evolution[END_REF] and the references therein).

1.3. Organization of the paper. In Section 2, we propose a general model for nonlocal heterogeneous diffusion, introducing in particular the notion of deciding factors, from which a wide range of nonlocal heterogeneous diffusion laws can be derived. Next, in Section 3, we draw a connection between the nonlocal and the local diffusion cases by a singular limit argument. This provides some new interpretation for equation (1.3). However, as was pointed out by two of the present authors in an earlier work [START_REF] Kim | Model for heterogeneous diffusion[END_REF], the diffusivity function D(x) may not be sufficient to describe heterogeneous diffusion accurately, because it does not fully cover velocity jump processes given by kinetic equations. In subsection 3.3, we make a related observation that the nonlocal diffusion equation, or position jump process, also leads to a range of local equations wider than (1.3). In Section 4, in order to answer the question "Where are the individuals?", we investigate steady states solutions (through both rigorous computations and numerical simulations) of the equations with heterogeneous nonlocal diffusion. This in particular reveals the role of deciding factors on the repartition of individuals, and may guide the use of these models in the applications. Finally, in Section 5, we briefly summarize our main contributions and present some perspectives.

Heterogeneous nonlocal diffusion equations

In a nonlocal ecological diffusion model, it is assumed that individuals have the ability to jump from one position to another at any given time.

In particular, we occasionally refer to the corresponding equation as the (position) jump process in order to distinguish it from kinetic equations which take integro-differential terms of velocity jumps. As mentioned earlier, a commonly used jump process is the homogeneous case (1.2), where the jump rate depends only on the distance between two points. However, if the environment is heterogeneous, the jump rate may depend on the departing point or arrival one. We may denote the jump rate by J(x, y) : the jump rate from x to y.

Let us emphasize again that this jump rate depends not only on the distance between the two points, but also on the direction of the jump. Thus the spatial heterogeneity of the environment along the way is counted by considering the distance and direction. This leads to a general nonlocal heterogeneous equation

(2.1) u t (t, x) = R N J(y, x)u(t, y) dy - R N J(x, y)u(t, x) dy.
Here, the first integral term stands for individuals arriving at x from all possible positions, whereas the second one accounts for individuals leaving from the position x. We refer to [START_REF] Hutson | The evolution of dispersal[END_REF] for its derivation from a discrete model in space and time. Notice that, provided that J is integrable with respect to each variable, the above equation is well-posed and satisfies the mass conservation property (at least formally) by Fubini's theorem.

Observe that (2.1) is recast

u t (t, x) = R N J(y, x)u(t, y) dy -m(x)u(t, x),
where (2.2) m(x) := R N J(x, y)dy stands for the total jump rate from position x. For late use, we also define

(2.3) J 1 (x) := R N
y -x J(x, y)dy the absolute first moment at position x. In other words, J 1 (x) stands for the average jump length among individuals jumping from the position x. Both the total jump rate and the average jump length are natural quantities which may be impacted by spatial heterogeneity. Obviously, equation (2.1) turns into the homogeneous nonlocal diffusion equation (1.2) by selecting J(y, x) = K(x -y) and defining m := R N K(z)dz. Yet our purpose is now to find practical ways to take into account the spatial heterogeneity and anisotropy in the choice of the jump rate function J. We will see that it is advantageous to specify J by distinguishing the effects of heterogeneity on the distance and on the direction of jumps.

2.1.

A few illuminating examples. Before we introduce general ways to handle heterogeneity, we start with some revealing examples. First, in the context of evolutionary biology, let us mention the work [START_REF]Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF] which considers a mathematical model for the fitness distribution in an asexual population under mutation and selection. Due to the presence of a phenotype optimum, the distribution of mutation effects on fitness depends on the parent's fitness (corresponding to x). Hence the authors need to consider "contextdependent" mutation kernels, corresponding in our setting to the case where (2.4) J(x, y) = K(x; y -x).

In other words a dispersal kernel K(x; •) has been assigned to each point x in the spatial domain, accounting for heteregeneity. In this above form, J(x, y) is the jump rate from x to y, and therefore the heterogeneity factor is considered at the departure point. This case arises from an Itô interpretation of a stochastic Poisson jump process, and therefore by analogy with the local diffusion we may call it nonlocal Itô type diffusion. To make the above form more meaningful, one may further specify the kernel function K, for instance by letting

(2.5) K(x; y -x) = m(x) K(y -x), or (2.6) K(x; y -x) = 1 g(x) N K y -x g(x)
,

where in either cases K denotes a probability density, which we may rescale without loss of generality so that its absolute first moment

R N z K(z)dz = 1.
In the former case, the x-dependence appears only in the factor m(x) which as outlined before stands for the total jump rate from the position x. In the latter case, one may compute that the total jump rate from x is constant equal to 1, while the average jump length (or absolute first moment at position x) is g(x) and now accounts for the heterogeneity. Both of these forms provide a tractable way to include heterogeneity in the model. However, depending on the context, one may want to consider a symmetrical situation where the jump rate from x to y is the same as the jump rate from x to y. The reasoning behind that may be that the jump rate of individuals may only depend on the path of the displacement, and not necessarily on heterogeneous conditions at the departure point. The above examples are ill-suited for such a situation, and leads us to instead consider the form

(2.7) J(x, y) = K x + y 2 ; y -x .
The symmetry of J(x, y) can then be reformulated as the evenness of the kernel K(p; •), for any p ∈ R N . Here p = x+y 2 stands for the point where heterogeneity dependence is involved: in this framework, heterogeneity is taken at the middle point of the jump.

The common trait to the above two forms (2.4) and (2.7) is that we assign a kernel to each position in the space domain. Then the heterogeneity corresponding to the path (x, y) is determined by environmental conditions at a single point p, which is then called a deciding factor. Extrapolating on the above examples, we are led to pick p = αx + βy, with α + β = 1, in other words as a barycenter of the departure and arrival points. When α = 1, we find the nonlocal Itô type diffusion where the influential environmental conditions are those at the starting point x. Let us notice that this particular case has been considered in [START_REF] Cortázar | A nonlocal inhomogeneous dispersal process[END_REF][START_REF] Cortázar | Existence and asymptotic behavior of solutions to some inhomogeneous nonlocal diffusion problems[END_REF][START_REF]Finite mass solutions for a nonlocal inhomogeneous dispersal equation[END_REF][START_REF]An inhomogeneous nonlocal diffusion problem with unbounded steps[END_REF]. When α = 1 2 we recover the symmetrical case when the jump rates are the same from both x to y and y to x. A key input of this work is to allow any barycenter of the starting and arrival points to be a deciding factor. This will be discussed again in subsection 2.3.

Still, it turns out that a single deciding factor may not be enough to fully understand the diffusion process, and we provide one last example to illustrate this. Let us consider, in spatial dimension N = 1, the case

(2.8) J(x, y) = 1 h(y) K y x 1 h(s) ds ,
where, as above, K is a probability density and h is positive. This nonlocal model arises from a Stratonovich (also referred to as Marcus [START_REF] Applebaum | Lévy Processes and Stochastic Calculus[END_REF]) integral of a Poisson jump process [START_REF] Bartlett | State-dependent jump processes: Itôstratonovich interpretations, potential, and transient solutions[END_REF]. Another way to derive it is to check that the heterogenous equation (2.1) with (2.8) reduces to the homogenenous nonlocal diffusion model (1.2) by to the change of variables

x = x 0 1 h(s)
ds.

This is consistent with the property that the Stratonovich integral preserves the chain rule of standard calculus. From the biological point of view, this means that the heterogeneity distorts the distances, as in the jump rate the (Euclidean) distance between x and y is replaced by the integral from x to y of the inverse of the function h. This integral may be understood as a distance with respect to a so-called food metric [START_REF] Choi | Chemotactic traveling waves by the metric of food[END_REF]. Due to the integral term in (2.8), one cannot isolate a single deciding point and it is the whole path which determines the heterogeneity. This partly explains why, in the next subsection, we will introduce a general framework where a kernel is assigned not only to any point (though, as suggested by the above examples, it is sufficient in some situations) but to any possible path. Such generality also allows us to consider the case of two deciding factors, which is convenient to distinguish the effects of heterogeneity on the total jump rate and on the average jump length.

2.2.

Assignment of a kernel to any possible path. In order to encompass all the above examples in a very general setting, our idea is the following: for any possible path, i.e. for any pair of points (x, y), we assign a dispersal kernel K(x, y; •). Notice that the case K(x, y; •) = K(x, •) corresponds to the situation in [START_REF]Dynamics of fitness distributions in the presence of a phenotypic optimum: an integro-differential approach[END_REF] mentioned above.

Precisely, we assume the form

J(x, y) = K(x, y; y -x), where K : (x, y; z) ∈ R N × R N × R N → K(x, y; z) ∈ R, is a nonnegative function such that, for any (x, y) ∈ R N × R N , m(x, y) := R N K(x, y; z) dz ∈ (0, +∞),
and

J 1 (x, y) := R N z K(x, y; z) dz ∈ (0, +∞).
Notice that, in particular, to any pair of points (x, y) there "corresponds" a nonlocal equation of type (1.2), where m(x, y) and K(x, y; •) are respectively the total jump rate and the dispersal kernel. Such a "fictive" equation could be interpreted as the diffusion law in a homogeneous medium where environmental conditions are everywhere the same as they are perceived by individuals jumping from x to y. In the sequel, even if "possibly fictive", we denote m(x, y) the total jump rate and J 1 (x, y) the absolute first moment, corresponding to the path (x, y).

With the above choice of the jump rate J, the general equation (2.1) is recast

(2.9) u t (t, x) = R N K(y, x; x -y)u(t, y) dy - R N K(x, y; y -x)u(t, x) dy.
To ensure that the right hand side of the above equation is "meaningful", we will also always assume that there exists

K ∈ L 1 (R N ) such that ∀(x, y) ∈ R N × R N , K(x, y; •) ≤ K.
Equations (2.1) and (2.9), together with the above assumptions, provide a general heterogeneous framework within which all the following sections will fit. Obviously, the latter equation (2.9) still includes the homogeneous equation: it naturally arises when we associate the same dispersal kernel to all paths, i.e. regardless of the departure and arrival points. It also includes (2.8), by letting

K(x, y; z) = 1 h(y) K z 1 0 1 h(x + s(y -x))
ds ,

provided that 0 < inf R N h ≤ sup R N h < +∞.
In order to emphasize the role of the total jump rate and the absolute first moment, we point out that a natural generalization of example (2.5) consists in

(2.10) K(x, y; z) = m(x, y) K(z),
where K is a probability density. In this case, the heterogeneity only appears in the total jump rate m(x, y) but not in the actual jump distribution K.

Similarly, a natural generalization of example (2.6) consists in

(2.11) K(x, y; z) = 1 g(x, y) N K z g(x, y) ,
where K is a probability density, rescaled so that its absolute first moment (2.12)

R N z K(z) dz = 1,
and the function g : R N × R N → R has positive infimum and supremum. Moreover, the total jump rate

m(x, y) = R N K(x, y; z) dz = R N K(z) dz
is actually spatially constant. However, the heterogeneity appears in the function g(x, y), which is the (fictive) relative absolute first moment of the dispersal kernel K(x, y; •), see (2.14) below, and can be interpreted as a (fictive) average jump length.

Finally, we may combine (2.10) and (2.11) so that both the total jump rate and the average jump length are spatially heterogeneous. The heterogeneous dispersal kernel can then be written as

(2.13) K(x, y; z) = m(x, y) g(x, y) N K z g(x, y) ,
where K is a probability density, rescaled so that its absolute first moment is equal to one. The relative absolute first moment of the dispersal kernel K(x, y; •) is defined by

(2.14) R N z K(x, y; z) dz R N K(x, y; z) dz
and is g(x, y) when assuming form (2.13).

2.3. The notion of deciding factors. In Section 3, we will compute the singular limits of (2.9) with focusing kernels to recover a wide range of local diffusion equations, including, but not limited to, (1.3). The point of such focusing kernels is that the average jump length goes to 0 as a singular limit parameter ε > 0 tends to 0. In this context, it turns out that it is mainly sufficient that the distribution kernel depends on a finite number of points. For instance, in the above example (2.13), one may assume that

g(x, y) = g(αx + βy), m(x, y) = m(α x + β y), where α + β = α + β = 1, so that J(x, y) = m(α x + β y) ( g(αx + βy)) N K y -x g(αx + βy) .
This means that the total jump rate m of the dispersal kernel corresponding to the path (x, y) is determined by environmental conditions at the single point α x + β y, which as outlined before we then call a deciding factor. On the other hand, the relative absolute first moment of the dispersal kernel corresponding to the path (x, y), which is also the average jump length g(x, y), is determined by environmental conditions at the (possibly) distinct deciding factor αx + βy. It is left to decide how to choose those two points, and to observe the consequences in the singular limit depending on that choice. For instance, one may pick α = 1 and β = 0, so that the average jump length g(x, y) is decided at the departure point x, while α = β = 1 2 , which means that the total jump rate m(x, y) is decided by the middle point x+y 2 . We believe that this interpretation, together with the connection between local and nonlocal diffusion which we next establish, provide a new perspective on the understanding of heterogeneous diffusion.

Singular limits: from nonlocal to local diffusion

In this section, we show how to obtain a large class of local diffusion laws from heterogeneous jump processes of the type (2.9). Recall that, in the homogeneous case, one can recover the heat equation (1.1) as the focusing kernel limit of the integro-differential (1.2). Precisely, assuming that K has a finite second moment the following holds: as ε → 0, the solution u ε of (1.2), with the kernel K(x) replaced by the focused kernel

K ε (x) := 1 ε N +2 K x ε ,
tends (in a sense we do not precise in this exploratory paper) to the solution of (1.1) with D depending on the second moment of K, starting with the same initial datum. This can be understood from a simple formal Taylor expansion that we will indicate below in a much more general setting. Hence, in some situations, nonlocal dipersal operators can be approximated by the local diffusion operator. This fact was long known: for instance, in evolutionary biology models, we refer to [START_REF] Kimura | A stochastic model concerning the maintenance of genetic variability in quantitative characters[END_REF] and [10, Chapter VI, subsection 6.4]. For a rigorous proof, one can use the explicit writing of the solution of (1.2) in the Fourier variable, see [START_REF] Andreu-Vaillo | Nonlocal diffusion problems[END_REF]Theorem 1.24]. As far as the singular limit of heterogeneous nonlocal diffusion equations is concerned, we refer to the recent works [START_REF] Molino | Nonlocal diffusion problems that approximate a parabolic equation with spatial dependence[END_REF][START_REF]Nonlocal approximations to Fokker-Planck equations[END_REF][START_REF] Sun | Approximate the Fokker-Planck equation by a class of nonlocal dispersal problems[END_REF] for some particular cases that all fall into our general framework below. Let us also notice the preprint [START_REF] Santos | A local/nonlocal diffusion model[END_REF] which starts from a local/nonlocal heterogeneous model.

In view of the very general framework outlined in Section 2, hereafter we will use the focusing kernels

K ε (x, y; z) := 1 ε N +2 K x, y; z ε .
Thanks to the "additional" variable z (ultimately replaced by y -x in the diffusion equation), the above form drives the average jump length to 0 as ε → 0, while preserving the scale of the heterogeneity in the singular limit.

3.1. The diffusivity matrix. The limiting local equation will involve the following notion of diffusivity.

Definition 3.1 (The diffusivity matrix). Let a dispersal kernel K(x, y; z) as in Section 2. For all p ∈ R N , we define a diffusivity matrix by

(3.1) D(p) := 1 2 R N z ⊗ zK(p, p; z) dz.
Alternatively, the diffusivity matrix D(p) = (d ij (p)) 1≤i,j≤n can be written component-wise as

(3.2) d ij (p) = 1 2 R N z i z j K(p, p; z) dz.
Notice that in the integral terms of the above definition, we only need to evaluate K(p, p ; z) at p = p . The reason is that the local equation arises in the limit of focusing kernels where the average jump length goes to 0, hence the departure and arrival point eventually coincide. There is however no mathematical obstacle to extend the above diffusivity matrix to the case where p = p . Notice also that it is implicitly assumed that the kernel K is such that

d ij (p) ∈ R for any p ∈ R N .
In order to illustrate this notion, let us consider the aforementioned example (2.13), namely K(x, y; z) = m(x, y) g(x, y) N K z g(x, y) ,

for which the diffusivity matrix is straightforwardly computed as

D(p) := m(p, p)g(p, p) 2 D.
Here, D is the constant diffusivity matrix associated with the homogeneous dispersal kernel K, while the total jump rate m and the average jump length g account for the environment's heterogeneities. Unsurprisingly, the diffusivity depends monotonically on both the total jump rate and the average jump length.

3.2.

The case of a single deciding factor. As outlined at the end of Section 2, we consider the nonlocal heterogeneous model (2.9) with

K(x, y; •) = K(αx + βy; •), where α + β = 1.
In other words, the nonlocal dispersal kernel associated with a path (x, y) is determined by the single point αx + βy. We immediately drop the tilde for convenience and rewrite (2.9) as (3.3)

u t (t, x) = R N K(αy + βx; x -y)u(t, y) dy -u(t, x) R N K(αx + βy; y -x) dy.
We point out that, when α = 0 and β = 1, the kernel is chosen from the arrival point and the equation is recast

u t (t, x) = R N K(x; x -y)u(t, y) dy -u(t, x) R N K(y; y -x) dy.
When α = 1 and β = 0, the kernel is chosen from the departing point and the equation is recast

u t (t, x) = R N K(y; x -y)u(t, y) dy -u(t, x)m(x),
where as before m denotes the integral of K with respect to its last variable. When α / ∈ {0, 1} the dispersal kernel is selected by a nontrivial linear combination of the departure and the arrival points. If we integrate (3.3) over x ∈ R N and use Fubini's theorem, we obtain as explained earlier the mass conservation property:

d dt R N u(t, x) dx = 0.
Here we implicitly assumed the framework introduced in Section 2: the function K(p; •) is nonnegative and bounded by some K ∈ L 1 (R N ) uniformly with respect to p ∈ R N . Moreover, we make the additional assumption that K is symmetric with respect to its last variable, i.e.

(3.4)

∀(p, z) ∈ R N × R N , K(p; -z) = K(p; z).
We also assume that K is "sufficiently smooth" in the p-variable, and that its derivatives with respect to the p-variable are "sufficiently z-integrable" (typically quadratic weights have to be supported). These assumptions play an important role in the singular limit argument, as will be clear from the terms appearing in the computations below.

Let us point out that a typical example of such dispersal kernels is

K(p; z) = m(p) K(z),
where K is a symmetric probability density, and thus only the mass m(p) of the dispersal kernel depends on position p. One may also assume that the relative absolute first moment, or average jump length, is spatially heterogeneous by considering

K(p; z) = m(p) g(p) N K z g(p)
.

We now formally derive local diffusion equations from the nonlocal diffusion equations (3.3). To do so, we actually take a sequence of focusing kernels. Denote, for any z ∈ R N , (3.5)

K ε (p; z) := 1 ε N +2 K p; z ε , 0 < ε 1.
Notice that these still satisfy our assumptions. However, with obvious notations, the total jump rate and the average jump length are related through

m ε (p) := R N K ε (p; z) dz = 1 ε 2 m(p), g ε (p) := R N z K ε (p; z) dz R N K ε (p; z) dz = εg(p).
Hence, as ε → 0, the average jump length goes to 0 but, in some sense, this is compensated by the fact that m ε (p) → +∞, i.e. individuals jump more often. In particular, the corresponding diffusivity matrix D ε of the focusing kernels is unchanged. This can be understood from the fact that, in the homogeneous case, the kernels (3.5) naturally arise from the change of variables

t ← ε 2 t , x ← εx.
Thus, in our context, focusing kernels can be understood as having the motion of individuals occur on a smaller scale than the heterogeneity. Let us proceed with our computation of the singular limit ε → 0. First, we plug the above focusing kernels into (3.3), use a change of variable z = y-x ε , and obtain

ε 2 u t (t, x) = R N K x + εαz; z u(t, x + εz)dz - R N K x + εβz; z u(t, x)dz,
where we used the symmetry of the dispersal kernel with respect to its last variable. Then

ε 2 u t (t, x) ≈ R N K(x; z) + εαz • ∇ p K(x; z) + ε 2 1 2 α 2 D 2 p K(x; z)z, z × u(t, x) + εz • ∇u(t, x) + ε 2 1 2 D 2 u(t, x)z, z dz - R N K(x; z) + εβz • ∇ p K(x; z) + ε 2 1 2 β 2 D 2 p K(x; z)z, z u(t, x)dz,
where ∇ p K(x; z) and D 2 p K(x; z) respectively denote the gradient vector and the Hessian matrix of the map of dispersal kernels, p → K(p; z).

Clearly the ε 0 order terms in the right hand side cancel each other. Moreover, since K(p; z) is symmetric with respect to z, we have

R N z • ∇ p K(x; z) dz = 0, R N K(x; z)z • v dz = 0,
for any v ∈ R N . Therefore, the ε 1 order terms in the right hand side also cancel out. As a result, recalling that β = 1 -α, we reach the local equation

u t (t, x) ≈ 1 2 R N K(x; z) D 2 u(t, x)z, z dz +α R N z • ∇ p K(x; z) z • ∇u(t, x) dz + 2α -1 2 R N D 2 p K(x; z)z, z u(t, x) dz.
We now define the "exponent" (3.6) q := 2 -2α (= 2β).

Then, using the exponent q and the coefficients d ij (p) of the diffusivity matrix from (3.2), we may rewrite the above diffusion equation as

u t (t, x) ≈ i,j d ij (x)∂ x i x j u(t, x) + (2 -q) i,j ∂ p j d ij (x)∂ x i u(t, x) +(1 -q) i,j ∂ p i p j d ij (x) u(t, x).
We pursue and, omitting to write the (t, x) variable, get

u t ≈ i,j ∂ x j d ij (x)∂ x i u +(1 -q) i,j ∂ p j d ij (x)∂ x i u + i,j ∂ p i p j d ij (x) u ≈ i,j ∂ x j d ij (x)∂ x i u + (1 -q) i,j ∂ x i ∂ p j d ij (x)u ≈ i ∂ x i j (d ij (x)∂ x j u + (1 -q)∂ p j d ij (x)u) ,
which is in the so-called divergence form. We finally obtain, in the singular limit, the (possibly anisotropic) diffusion equation

(3.7) u t = i ∂ x i j d q ij (x)∂ x j d 1-q ij (x)u ,
which is an anisotropic generalization of the form (1.3).

Remark 3.2 (Isotropic diffusion)

. The relation between (3.7) and (1.3) becomes obvious in the isotropic case, i.e. when radial symmetry is assumed for the mapping z → K(p; z), which is a stronger assumption than (3.4). In this case, we obtain D(p) = D(p)I, where I is the identity matrix and

D(p) := 1 2N R N z 2 K(p; z) dz.
The diffusion equation (3.7) is then recast

u t = ∇ • D q (x)∇ D 1-q (x)u , which is exactly (1.3).
Let us first observe that, if q = 1 in (3.7), we reach Fick's law:

u t = ∇ • (D(x)∇u). (q = 1, α = 1
2 ) This is the case when α = β = 1 2 , meaning that the middle point between departure and arrival points is used to define the spatial heterogeneity in the nonlocal model (for instance, in a heterogeneous total jump rate m).

The next case is when q = 0 in (3.7), which yields Chapman's law:

u t = ∇ • (∇ • (D(x)u)). (q = 0, α = 1)
Here, ∇ • (Du) is the vector that consists of the divergence of the rows of the matrix Du, that is

∇ • (Du) = ∇ • (d i1 u, ..., d iN u) 1≤i≤N .
This is the case when α = 1, β = 0, meaning that only the departure point is used to define the spatial heterogeneity in the nonlocal model. For a general q, we may rewrite equation (3.7) as follows. First, observe that

j ∂ x j d ij u = j d q ij ∂ x j d 1-q ij u + j ∂ x j (d q ij )d 1-q ij u.
Let N q (x) = ((N q ) i (x)) be the vector given by

(3.8) (N q ) i (x) = j ∂ x j (d q ij (x))d 1-q ij (x) .
Then, the component-wise equation (3.7) can be written as (3.9)

u t = ∇ • ∇ • (D(x)u) -N q (x)u .

Remark 3.3 (Orthotropic diffusion).

Assume the component-wise symmetry, which is a stronger than the symmetry assumption (3.4): for all

(p, z) ∈ R N × R N , (3.10) K(p; α 1 z 1 , ..., α N z N ) = K(p; z) for all (α i ) ∈ {±1} N .
Then, the diffusivity matrix D(p) becomes diagonal, which is often called orthotropic. Under assumption (3.10), equation (3.7) turns into

u t = i ∂ x i d q ii (x)∂ x i d 1-q ii (x)u .
Furthermore, since the diffusivity matrix D(p) is diagonal, D q (p) can be defined element-wise, and the equation is recast where

u t = ∇ • D q (x) ∇ • (D 1-q (x)u) .
α + β = α + β = 1,
and ν, K will be precised below. As explained at the end of Section 2, this means that the heterogeneity in the nonlocal dispersal kernel associated with a path (x, y) is no longer determined by a single point, but by the two points αx + βy and α x + β y. As in the previous subsection, both points are still chosen as linear combinations of the departure and arrival points. Notice that the diffusivity matrix associated with a dispersal kernel of the form (3.11) remains the same, regardless of whether α = α or not, that is regardless of whether the two deciding factors are actually the same or not. Therefore one may expect that the limit local equation is also unchanged. As we will see below, it turns out that this is not the case. In some sense, the diffusivity matrix does not convey all the information of the motion in the heterogeneous case. This is quite similar to the observation made by two of the present authors in [START_REF] Kim | Model for heterogeneous diffusion[END_REF] in the kinetic model.

The motivation behind such a double heterogeneity is the fact that the dispersal kernel K(x, y; •) can be mostly characterized by its mass (or zero moment) m(x, y) := R N K(x, y; z) dz, and its relative absolute first moment

g(x, y) := R N z K(x, y; z) dz R N K(x, y; z) dz .
Indeed, as we explained in Section 2, the former can be interpreted as the total jump rate of individuals, and the latter as the average jump length. Both notions are crucial in the modelling of diffusion in ecology. It is then natural to allow these two functions not only to depend on space, but also to vary independently and be decided at possibly distinct points.

A typical example of dispersal kernel in the form (3.11) is provided by

(3.12) K(x, y; z) = m(α x + β y) g(αx + βy) N K z g(αx + βy) ,
where K is a fixed probability density, rescaled so that (2.12) holds.

We now proceed with the general form (3.11). Dropping again the tilde for convenience, the nonlocal equation (2.9) rewrites (3.13)

u t (t, x) = R N ν(α y + β x)K(αy + βx; x -y)u(t, y) dy -u(t, x) R N ν(α x + β y)K(αx + βy; y -x) dy.
According to (3.11), it is assumed that the two variables accounting for the heterogeneity are separable in the dispersal kernel. Moreover, in order to satisfy the assumptions of Section 2, we take ν : R N → R to have positive infimum and supremum, and K : R N ×R N → R a nonnegative and nontrivial function. We also assume that K(p; z) ≤ K(z) for some K ∈ L 1 (R N ). As in the previous section, we impose the symmetry condition,

∀(p, z) ∈ R N × R N , K(p; -z) = K(p; z),
that both ν(p ) and K(p; z) are "sufficiently smooth", and that all derivatives of ν(p )K(p; •) in the p, p -variables are "sufficiently z-integrable".

We now consider the focusing kernels

(3.14) ν ε (p ) := 1 ε 2 ν(p ) K ε (p; z) := 1 ε N K p; z ε , 0 < ε 1.
Plugging these into the nonlocal equation, and by a change of variable z = y-x ε , we obtain

ε 2 u t (t, x) = R N ν(x + εα z)K(x + εαz; z)u(t, x + εz) dz - R N ν(x + εβ z)K(x + εβz; z)u(t, x) dz.
After performing formal asymptotic expansions and using similar arguments as in subsection 3.2, we reach

u t (t, x) ≈ ν(x) R N K(x; z) 1 2 D 2 u(t, x)z, z dz +αν(x) R N z • ∇ p K(x; z) z • ∇u(t, x) dz +α R N (z • ∇ν(x))K(x; z)(z • ∇u(t, x)) dz + (2α -1)ν(x) R N 1 2 D 2 p K(x; z)z, z dz +(α + α -1) R N (z • ∇ν(x))(z • ∇ p K(x; z)) dz +(2α -1) R N 1 2 D 2 ν(x)z, z K(x; z) dz u(t, x). Denote q := 2 -2α = 2β, q := 2 -2α = 2β .
In the setting (3.11) and after dropping the tilde, we observe that the coefficients d ij (p) of the diffusivity matrix D(p) of Definition 3.1 are given by

d ij (p) = ν(p) 1 2 R N z i z j K(p; z) dz.
In light of this, proceeding as in subsection 3.2, we can rewrite the above equation in a conservative form as

(3.15) u t = i ∂ x i j ν q -q (x)d q ij (x)∂ x j ν q-q (x)d 1-q ij (x)u .
Notice that there are two degrees of freedom, given by q and q . First, find that, if q = q , then (3.15) is written as

(3.16) u t = i ∂ x i j d q ij (x)∂ x j d 1-q ij (x)u .
In this case, diffusivity in the sense of Definition 3.1 alone decides the diffusion phenomenon. This is exactly the situation considered in subsection 3.2 with "a single deciding factor". For instance, we reach

u t = ∇ • (∇ • (D(x)u)), (q, q
) = (0, 0)

u t = ∇ • (D(x)∇u).
(q, q ) = (1, 1).

However, as seen in (3.15), this is no longer the case when two different types of heterogeneity are considered.

Example 3.4. Let us return to our previous example (3.12) to better understand the meaning and the role of ν: when

K(x, y; z) = m(α x + β y) g(αx + βy) N K z g(αx + βy) ,
then ν = m is the total jump rate. Therefore the limiting local equation (3.15) typically involves not only the diffusivity matrix but also the total jump rate of the original nonlocal model. Let us also recall that, in this example, the diffusivity matrix is

D(p) = m(p)g(p) 2 D, where D = ( d ij ) ij is spatially homogeneous. Then (3.15) rewrites (3.17) u t = i ∂ x i m q (x)g 2q (x) j d ij ∂ x j m 1-q (x)g 2-2q (x)u .
In this form, one may identify the role of not only the total jump rate m but also the average jump length g. In particular, the parameters q and (α, β) correspond to the average jump length g while q and (α , β ) correspond to the total jump rate m.

Example 3.5. Let us first consider the case when α = 1 and α = 1/2.

In the above example, this means that the total jump rate is determined by the departure point, while the average jump length is determined by local conditions on the travel path, which can be reasonably simplified by choosing the middle point. From (3.15), we get

u t = ∇ • (ν -1 (x)D(x)∇(ν(x)u)). (q, q ) = (1, 0)
Yet, we may instead consider the opposite case where the average jump length is determined at the departure point while the total jump rate is determined at the middle point (for instance, individuals may choose to go back to the departure point if they see a decline in the environmental conditions). This corresponds to α = 1/2 and α = 1. From (3.15), we get

u t = ∇ • (ν(x)∇ • (ν -1 (x)D(x)u)), (q, q ) = (0, 1)
This equation is similar to a local equation derived in [START_REF] Kim | Model for heterogeneous diffusion[END_REF] from a velocity jump process. The only difference is that ν here should be µ -1 in [START_REF] Kim | Model for heterogeneous diffusion[END_REF]. To make it exactly the same, one needs to choose q = 0, q = -1 so that we get

u t = ∇ • (ν -1 (x)∇ • (ν(x)D(x)u)). (q, q ) = (0, -1)
Still, from the point of view of the position jump process, the case (q, q ) = (0, 1) seems more reasonable. Indeed, it is clear that diffusivity here is increasing with respect to ν, which is typically the total jump rate. However, in the kinetic model of [START_REF] Kim | Model for heterogeneous diffusion[END_REF] the function µ stands for the turning frequency, with respect to which the diffusivity was decreasing. Thus, both notions should not be confused and it is not surprising that positive and negative exponents are swapped.

For general q and q , there are many ways to write the limiting local equation (3.15). We decide to write it as a "perturbation" of the case (q, q ) = (0, 1) considered in the above example. To do so we write

j ν∂ x j ν -1 d ij u = j ν q -q d q ij ∂ x j (ν q-q d 1-q ij u) + j ∂ x j (ν q -q-1 d q ij )ν 1+q-q d 1-q ij u = j ν q -q d q ij ∂ x j (ν q-q d 1-q ij u) +νu(N q,q ) i ,
where N q,q (x) = ((N q,q ) i (x)) is a correction vector given by

(3.18) (N q,q ) i := j ∂ x j (ν q -q-1 d q ij )ν q-q d 1-q ij ,
which vanishes in the case (q, q ) = (0, 1). As a result, the diffusion equation (3.15) is recast

(3.19) u t = ∇ • ν(x)∇ • (ν -1 (x)D(x)u) -ν(x)N q,q (x)u .
u t = i ∂ x i ν q -q (x)d q ii (x)∂ x i ν q-q (x)d 1-q ii (x)u .
Furthermore, since the diffusivity matrix D(p) is diagonal, D q (p) can be defined element-wise. Then, we may write it as

(3.20) u t = ∇ • ν q -q (x)D q (x) ∇ • (ν q-q (x)D 1-q (x)u) . Remark 3.7 (Isotropic diffusion). If z → K(p; z) is radial symmetric, we obtain
D(p) = D(p)I, D(p) := ν(p) 2N R N z 2 K(p; z) dz.
The isotropic diffusion equation is now

(3.21) u t = ∇ • ν q -q (x)D q (x)∇ ν q-q (x)D 1-q (x)u .
For instance, if (q, q ) = ( 1 2 , 1), we obtain

u t = ∇ • ν(x)D(x) ∇ ν -1 (x)D(x) u .

Steady states and numerical simulations

In this section we present some enlightening observations regarding the steady states of (2.1), when the jump rate is of the form This connects with the issue discussed in, e.g., [START_REF] Cantrell | Evolutionary stability of ideal free dispersal strategies: a nonlocal dispersal model[END_REF][START_REF] Cosner | Evolutionary stability of ideal free nonlocal dispersal[END_REF] and is complemented by some numerical simulations. Those simulations are performed on a bounded domain, B(0, R) ⊂ R N , see subsection 4.2 for details.

4.1. Steady states. We consider the case of a single deciding factor with J(x, y) of the form (4.1). In what follows, we omit the tilde notation for readability. The nonlocal model (2.1) is then recast (4.2)

u t = R N [m(αy + βx)u(t, y) -m(αx + βy)u(t, x)] K(y -x)dy, α + β = 1.
As seen in subsection 3.2, taking the focusing kernel limit leads, in the isotropic case, to the local model

(4.3) u t = ∇ • D q (x)∇ D 1-q (x)u , q = 2 -2α, with D(x) = 1 2 m(x)k, where k := 1 N R N z 2 K(z)dz.
For the rest of this section, we consider a nonnegative, nontrivial initial data u 0 ∈ L 1 (R N ). Also, we assume that m is a nonnegative continuous function. Finally, for numerical simulations and unless otherwise stated, we shall use the Gaussian kernel

K(z) := C N e -a N z 2 ,
where C N > 0, a N > 0 are such that not only R N K(z)dz = 1 but also R N z K(z)dz = 1. In particular C 1 = a 1 = 1 π . The case α = 0 (q = 2). In this case, we observe that u(x) ≡ Cm(x) is a steady state of (4.2) for any constant C ≥ 0. This is consistent with the local counterpart, that is (4.3), for which u(x) ≡ CD(x) is also a steady state. In other words, at equilibrium, the population is mainly localized at positions x where m(x), or equivalently D(x), is large. This is in agreement with the model since the total jump rate from x to y is by m(y), meaning the jump rate is higher for jumps that arrive at high values of m, and conversely smaller for jumps that arrive at low values of m. Therefore the solution of (4.2) is expected to converge, at large time, to

(4.4) p(x) := Cm(x), C = R N u 0 (x)dx R N m(x)dx ≥ 0,
the value of C ensuring the mass conservation. In particular, if m / ∈ L 1 (R N ), we expect that p(x) ≡ 0. We now present some numerical simulations to illustrate the above. Since these simulations are done on the bounded domain Ω R = B(0, R) ⊂ R N , the expected profile is (see subsection 4.2 for details)

(4.5) p R (x) := C R m(x), C R = Ω R u 0 (x)dx Ω R m(x)dx > 0,
which is consistent with our simulations, see Figure 1. Moreover, we readily check that p R → p uniformly on compact sets of R N , and uniformly on R N for bounded m. However, for general m, this convergence may not be uniform on R N , as can be seen with N = 1, m(x) = e x , for which ||p R || ∞ → R u 0 (x)dx while p ≡ 0.

The case α = 1 (q = 0). Here, we assume m(x) > 0 for all x. We observe that u(x) ≡ Cm(x) -1 is a steady state of (4.2) for any constant C ≥ 0. This is again consistent with the local counterpart, that is (4.3), for which u ≡ CD(x) -1 is also a steady state. The situation is thus reversed compared to the case α = 0: the population will mainly concentrate at positions x where m(x) is small. This is, again, in agreement with the model since the jump rate from x to y is given by m(x), meaning the jump rate is higher for jumps that depart from high values of m, and conversely smaller for jumps that depart from low values of m. Therefore the solution of (4.2) is expected to converge, at large time, to

(4.6) p(x) := Cm(x) -1 , C = R N u 0 (x)dx R N 1 m(x) dx ≥ 0,
the value of C ensuring the mass conservation. In particular, if m -1 / ∈ L 1 (R N ), we expect that p(x) ≡ 0. Similarly to the case α = 0, on the bounded domain Ω R = B(0, R) ⊂ R N , the expected profile would be

(4.7) p R (x) := C R m(x) -1 , C R = Ω R u 0 (x)dx Ω R 1 m(x) dx > 0,
which is consistent with our simulations, see Figure 2. As in the case α = 0, we readily check that p R → p uniformly on compact sets of R N .

The case α = 1 2 (q = 1). In that setting, it is clear that constants are steady states of both (4.2) and (4.3). Therefore the profile of m plays no role in the distribution of population at equilibrium. This is again expected since then

J(x, y) = m x + y 2 K(y -x) = J(y, x),
that is the jump rate from x to y is equal to that from y to x. Therefore the solution of (4.2) is expected to converge, at large time, to p(x) ≡ 0, due to the mass conservation. Note that on a bounded domain Ω R = B(0, R) ⊂ R N , we however expect a convergence towards the constant p R = 1 |Ω R | Ω R u 0 (x)dx, which is in agreement with our numerics, see Figure 3. Clearly, p R → 0 as R → +∞.

Towards general 0 ≤ α ≤ 1. The steady states are thus well-understood for general functions m in the cases α ∈ {0, 1 2 , 1}. In particular, those steady states depend only on the total jump rate m, and not on the jump distribution K. We aim at extending those results for general α. Hence, for a given α ∈ {0, 1 2 , 1} and a given nonnegative and nontrivial m, we look for a nonnegative and nontrivial u such

(4.8) m(αy + βx)u(y) = m(αx + βy)u(x), ∀x, y ∈ R N ,
which makes u = u(x) a steady state for (4.2). Notice that such a function u is not ensured to exist, as can be seen with m(•) = || • ||: in this case, for x = 0, (4.8) enforces u(y) ≡ β α u(0) for any y = 0; returning to (4.8), this is absurd if α = 1 2 . However, there are some profiles m for which such a function u exists for any 0 ≤ α ≤ 1. Indeed, assuming N = 1, m(x) = e ax with a ∈ R, we see that for each α, the function u α (x) = e (β-α)ax = m(x) β-α solves (4.8). Because u α / ∈ L 1 (R), we expect the solution of (4.2) to converge towards p α (x) ≡ 0 due to the mass conservation. However, on the bounded domain . In solid lines, for various α, the numerical approximation of u at t = 4000, which all coincide with p α R .

Ω R = [-R, R], the expected profile is (4.9) p α R (x) := C α R u α (x), C α R = Ω R u 0 (x)dx Ω R u α (x)dx > 0,
which is consistent with our simulations, see Figure 4, for which we have chosen u 0 ≡ m which, obviously, is not in L 1 (R). Nevertheless, our goal here is to illustrate the fact that, as α increases, there is a smooth transition from m (when α = 0) to m -1 (when α = 1), the "switch" precisely occurring at α = 1 2 . Last, one can extend this example to N ≥ 2 with m(x) = e a,x and a ∈ R N , for which u(x) = e (β-α) a,x solves (4.8).

Another example with N = 1 is given by the Gaussian profile m(x) = e -ax 2 with a > 0, for which u α (x) = e -(β 2 -α 2 )ax 2 solves (4.8). Since u α ∈ L 1 (R) if and only if α < 1 2 , we expect that the solution of (4.2) converges towards

p α (x) = C α u α (x) if α < 1 2 , 0 if α ≥ 1 2 , C α = R N u 0 (x)dx R N u α (x)dx > 0. Meanwhile, on the bounded domain Ω R = [-R, R], since u α ∈ L 1 (Ω R
) for all α, we expect that the solution converges towards p α R given by (4.9). This is what we observe numerically, see Figure 5. Note that the above comments for the exponential profile still apply here. Last, one can extend this example to N ≥ 2 with m(x) = e -a,x 2 , where a ∈ (0, +∞) N and x 2 := (x 2 1 , . . . , x 2 N ), for which u(x) = e -(β 2 -α 2 ) a,x 2 solves (4.8).

Where are the individuals? The role of α. Our above observations and results can be summarized as follows: when α = 0 (resp. α = 1, resp. α = 1 2 ), some steady states are proportional to m (resp. m -1 , resp. x → 1). We 100 . In solid lines, for various α, the numerical approximation of u at t = 4000, which all coincide with p α R . therefore conjecture that, as α increases, the population should concentrate at positions x where m(x) is small. This was proved above for the particular cases of exponential and Gaussian profiles of m. To further support this conjecture, we investigate the case of a two-patch asymmetric function m, namely (4.10) m(x) = γ(x + 5) + 2γ(x -5), γ(x) = 1 1 + x 2 . The numerical solutions are displayed in Figure 6, for different values of α. We observe indeed that, as α increases, the population shifts away from the positions where m is large and settles in positions where m is small. The role of the tails of K. In all the above simulations, we always considered that the probability density K (rescaled so that R N ||z||K(z)dz = 1) has Gaussian tail, that is (say N = 1)

K(z) = 1 π e -1 π z 2 =: K 1 (z).
In order to draw a comparison, we also consider the following kernels: In solid lines, for various α, the numerical approximation of u at t = 4000 . The initial data is of mass Ω R u 0 (x)dx = 4.

hold and that K 4 is not radially decreasing. In particular, for i ≥ 2, all the kernels K i have the same second moment, which means that they all lead to the same local model in the focusing kernel limit.

Our goal here is to numerically compare the steady states depending on the chosen kernel. When α ∈ {0, 1 2 , 1}, we do not observe, as expected, any difference, for various functions m. The situation is quite different when, say, α = 1 4 , see Figure 7. Indeed, we observe the following: the heavier the tails of K i , the more the corresponding solution u i tends to concentrate at positions where m is large, as can be seen from the first and last zoomed parts. Knowing that, one may expect that u j -u i changes sign only once on (0, +∞). It turns out to be the case for (i, j) = (4, 3) but not for (i, j) = (2, 3), see Figure 8. The latter also indicates that, in a neighbourhood of x = 0, u 3 -u 2 is (positive and) very small while u 4 and u 3 are clearly distinct.

4.2. Numerical implementation. Since any numerical computation has to be done on a bounded domain Ω ⊂ R N , the numerical integration of the right-hand side of (2.1) is not trivial, even for simple cases such as J(x, y) = K(y -x). Indeed, in order to compute the convolution K * u close to the boundary ∂Ω, one would have to extrapolate u outside of Ω itself. To circumvent this, we consider the following problem on the bounded domain Ω ⊂ R N : In other words, we assume that any jump that would leave Ω is in fact omitted. The functions p R and p α R appearing in subsection 4.1 are in fact steady states of (4.11) with Ω = Ω R = B(0, R) ⊂ R N , and we have checked Here u i corresponds to the numerical solution when K = K i , in the conditions of Figure 7.

that our numerical solution converge towards them for some large values of R. As mentioned in subsection 4.1, the steady states p R and p α R typically converge, as R → +∞, to a steady state of (4.2), locally uniformly in R N .

Summary and perspectives

For the modelling of multi-dimensional nonlocal heterogeneous diffusion, we have proceeded as follows: to any possible path (x, y) we have assigned a dispersal kernel K(x, y; •) such that K(x, y; -z) = K(x, y; z), and relevant quantities for this path are m(x, y) := We have revealed, in the focusing kernel limit, the connection Nonlocal model (3.3) q=2-2α

-----→ diffusion equation (3.7), or (3.9),

3. 3 .

 3 The case of two deciding factors. We now consider the nonlocal heterogeneous model (2.9) with(3.11) K(x, y; •) = ν(α x + β y) K(αx + βy; •),

Remark 3 . 6 (

 36 Orthotropic diffusion). Under the assumption (3.10) for orthotropic diffusion, equation (3.15) turns into

(4. 1 )

 1 J(x, y) = m(αx + βy) K(y -x).

Figure 1 . 1 100 x 2 .

 112 Figure 1. Here α = 0 and R = 10. In dotted black, the function m. In blue, the numerical approximation of u at time t = 4000. In dotted red, the expected profile p R (x). On the left m(x) = 1 1+x 2 , and on the right m(x) = 1 + 1 100 x 2 . The outcome depends only on the initial mass, here Ω R u 0 (x)dx = 4.

Figure 2 .

 2 Figure 2. Conditions are the same as in Figure 1, except α = 1.

Figure 3 .

 3 Figure 3. Conditions are the same as in Figure 1, except α = 1 2 .

Figure 4 .

 4 Figure 4. In dotted black, the function u 0 (x) = m(x) = e ax with a = ln 2 10. In solid lines, for various α, the numerical approximation of u at t = 4000, which all coincide with p α R .

Figure 5 .

 5 Figure 5. In dotted black, the function u 0 (x) = m(x) = e -ax 2 with a = ln 100100 . In solid lines, for various α, the numerical approximation of u at t = 4000, which all coincide with p α R .

K 2 ( 4 .

 24 Note that, denotingk i = R z 2 K(z)dz, we have k 1 = π 2 and k 2 = k 3 = 2.Finally, we consider one last kernel K 4 with a heavier tail, namelyK 4 (z) := 1 1 + |z| 3 A log(2 + |z|) 3 -B log(2 + |z|) 5 + C log(2 + |z|) 7 ,where A ≈ 2.32, B ≈ 2.38 and C ≈ 0.62 are numerically selected such that R K 4 (z)dz = 1, R |z|K 4 (z)dz = 1 and k 4 = 2. Notice that K 4 ≥ 0 does

Figure 6 .

 6 Figure 6. In dotted black, the function m given by (4.10).In solid lines, for various α, the numerical approximation of u at t = 4000 . The initial data is of mass Ω R u 0 (x)dx = 4.

  αy + βx)u(t, y) -m(αx + βy)u(t, x)] K(y -x)dy.

Figure 7 . 1 1+x 2 .

 712 Figure 7. Here α = 1 4 . In dotted black, the function m(x) = 1 1+x 2 . In red, the numerical approximation of the solution at t = 10 4 , corresponding to the kernel K 1 (red), K 2 (yellow), K 3 (blue) and K 4 (purple). Some parts of the first image have been zoomed in to improve visibility. On the first and second zooms, the yellow curve almost coincides with the blue one. The initial data is of mass Ω R u 0 (x)dx = 4. The results are similar on larger domains.

Figure 8 .

 8 Figure 8. Left: u 4 -u 3 . Right: u 3 -u 2 .Here u i corresponds to the numerical solution when K = K i , in the conditions of Figure7.

  , y; z) dz the total jump rate, and g(x, y) := R N z K(x, y; z) dz R N K(x, y; z) dz the average jump length. Also, for p ∈ R N , we have defined D(p) := d ij (p) := 1 2 R N z i z j K(p, p; z) dz 1≤i,j≤n the diffusivity matrix. 5.1. A single deciding factor. First we have considered the case K(x, y; z) ← K(αx + βy; z), α + β = 1.

which is possibly anisotropic. In the so-called orthotropic case, the diffusivity matrix is diagonal and the limit diffusion equation is recast

In the so-called isotropic case, the diffusivity matrix is scalar and the limit diffusion equation is recast

The following examples recover some standard diffusion laws and are worth being mentioned.

Deciding factor

Limit diffusion equation Limit law

Furthermore, in the framework K(αx+βy; z) ← m(αx+βy)K(z), we have investigated the form of the steady state solutions. This revealed that, as α increases, the population shifts away from the positions where m is large and settles in positions where m is small. In some situations, the steady states are fully determined by α and m. Nevertheless, when α / ∈ {0, 1 2 , 1}, we have tested the influence of different kernels K and noticed some subtle differences on the shape of the steady states. It seems that, the heavier the tails of K, the more the solution concentrates, which is slightly counter-intuitive. Understanding how the steady states are determined by the interplay of α, m and K appears very challenging and deserves further investigations.

Two deciding factors. Next we have considered the case

We have revealed, in the focusing kernel limit, the connection Nonlocal model (3.13)

which is possibly anisotropic. In the so-called orthotropic case, we refer to Remark 3.6. In the so-called isotropic case, we refer to Remark 3.7 and the following examples are worth being mentioned (where D.F. means deciding factor). D.F. for average jump length D.F. for total jump rate Limit diffusion equation

In particular this reveals that, in the case of two distinct deciding factors, the diffusivity matrix is not enough to characterize heterogeneous diffusion. More precisely, the choice of two deciding factors allows us to recover the whole range of local diffusion equations that has been derived from kinetic equations in [START_REF] Kim | Model for heterogeneous diffusion[END_REF]. 5.3. Many or infinitely many deciding factors. The two above subsections suggest that, as the number of deciding factors increase to account for more general heterogeneity, one may recover new local models through the singular limit procedure. In particular, it is worth pointing out that Wereide's law does not only appear as a limit of a single deciding factor "close to departure" (see subsection 5.1), but also as the singular limit of the focusing kernels associated with (2.8). We did not detail this computation since it proceeds similarly to what we have presented. In other words, the nonlocal Stratonovich model does converge to the Stratonovich type diffusion, as one may expect. We believe that generalizations of (2.8), for which the whole path decides, deserve further investigations.


Mohammed El

Amin Tebib

Pascal André
email: pascal.andre@univ-nantes.fr

Olivier Cardin
email: olivier.cardin@univ-nantes.fr

André Pascal

Cardin Olivier

A Model Driven Approach for automated generation of Service-oriented Holonic Manufacturing Systems

Keywords: Service-Oriented Architecture, Holonic Manufacturing Systems, Model-Driven Engineering, Verification

In the context of manufacturing in Industry 4.0, software systems become of prime importance. Efficient, adaptable and trusted software services are required. Several approaches succeeded in creating a Service-oriented Holonic Manufacturing System that combines the advantages of Service-oriented Architectures and Holonic Manufacturing Systems. However these systems until now suffer from many shortcomings, among which genericity, lack of proof of the functional behaviour correctness, architecture modularity, etc. These systems are often manually implemented and become hardly adaptable and reconfigurable to different contexts (resources, workshop...). In this paper, we investigate a Model Driven Engineering approach to represent these systems in order to automate the generation of the services logic code from an abstract models and construct a new software chain that deals with all the shortcomings cited above.

Introduction

Industry 4.0 increases the influence of software in the field of cyber-physical systems (CPS) where the service paradigm becomes pregnant. Cyber-Physical Production Systems (CPPS) relies on the latest and foreseeable further developments of computer science (CS), information and communication technologies (ICT), and manufacturing science and technology (MST). CPPS consist of autonomous and cooperative elements and subsystems that are connected based on the context within and across all levels of production, from processes through machines up to production and logistics networks [START_REF] Monostori | Cyber-physical systems in manufacturing[END_REF]. Business processes can be connected to manufacturing processes. Cyber manufacturing is a transformative concept that involves the translation of data from interconnected systems into predictive and prescriptive operations to achieve resilient performance [START_REF] Lee | Introduction to cyber manufacturing[END_REF].

The unifying paradigm between processes and systems is the concept of service, including the cloud stack XaaS and the Service Oriented Architectures (SOA). The service paradigms supports scalability from (high level) business processes in enterprise architecture frameworks to (low level) hardware operations. The coordination and control of such complex systems by the way of actors or components requires methods and techniques to design, verify and deploy the services, possibly on the fly, making service engineering an unavoidable approach to develop new generation cyber-manufacturing systems. This is part of requirements of CPS mentioned by Wang et al. in the category 'design methodology' [START_REF] Wang | Current status and advancement of cyber-physical systems in manufacturing[END_REF]. As mentioned by Bauer et al., the traditional automation pyramid is dissolving and manufacturing IT is moving towards service-orientation and app-orientation [START_REF] Bauer | Movement towards serviceorientation and app-orientation in manufacturing {IT}[END_REF]. As an example, in cloud manufacturing, SoA was identified to meet the requirements of all higher level manufacturing CPS layers due to the reduced time constraints present [START_REF] Morgan | The cyber physical implementation of cloud manufactuirng monitoring systems[END_REF].

However "servicizing" is not sufficient to get software product with good-enough quality as defined by Meyer [Meyer, 1988]. First, we need software development techniques to improve both the quality of the resulting applications (correctness and robustness) but also its maintenance (extendibility, reusability). Second, we need analysis tools to check the service model properties on various aspects (structure, dynamics, functional and non-functional), and model transformations to compute new models or to generate code in the spirit of Model Driven Engineering (MDE).

In a previous article [START_REF] André | Trusted Services for Cyber Manufacturing Systems[END_REF], we advocated for a new vision of CPPSs construction with an emphasis on trusted service based component systems to represent the software part of CPPSs. The main idea was to incorporate a broad set of software design techniques and practices to improve the quality of the delivered CPPSs. In this article, we follow this position when revisiting a legacy SOA oriented CPPSs application and showing by practice our vision. We mainly address the software construction process and maintenance issues (modularity, extendibility, reconfigurability, evolution) than the correctness and robustness issues which details will be postponed to future contributions. Legacy applications miss abstraction, genericity and modularity which make harder the evolution or reuse of such applications. Thanks to model driven engineering and software development best practice, we propose a software construction process that improves abstraction and separates the concerns to face the above drawbacks. Abstraction enables to reason at the good level e.g. verification of model properties before implementation. The separation of concerns enables to separate the problem domain from the implementation issues (simulation, real workshop, user interface, interoperability and communications...)

We present the problem statement in section 2 and sketch an overview of our proposal in section 3. Thanks to MDE, building the manufacturing software becomes itself a Software Product Line where we intend to rationalize the software production and maintenance, including features and variability, shared aspects are separated from specific aspects. The approach is illustrated by on-going experimentation on a benchmark workshop in section 4. The application of this work is discussed in section 5. In conclusion, we draw the open perspectives and expected advantages.

Problem Statement

The starting point is a double finding, from literature review and current practice.

Related works mention that service engineering in the context of CPPSs is still a craft activity, usually at the implementation level [START_REF] Rodrigues | Self-interested service-oriented agents based on trust and qos for dynamic reconfiguration[END_REF], Morariu et al., 2013]. Two main levers are still required to go further. First, we need service models that can fit to various semantics and various granularity levels. Indeed, the concept of CPPS covers many classes of (physical) systems, from the manufacturing workshop to the whole supply chain. Taking the example of HMS, the control of such systems is often recursive, if not fractal, in order to aggregate the available resources and enable a heterarchic control architecture. Therefore, services that might be used at various levels of the architecture need to fit various granularity and the portability of services between different applications with their own semantic requires an adaptability of the services to be effective.

In practice, we started from the work of [START_REF] Quintanilla | A modeling framework for manufacturing services in service-oriented holonic manufacturing systems[END_REF] who proposed a new vision for the conception of agile HMS by coupling HMS architectures with Service-oriented ones, which gave birth to the concept of Service-oriented Holonic Manufacturing Systems (SoHMS) [Quintanilla, 2015]. They implemented a manufacturing software in Java to control an assembly line as illustrated by Figure 1. The Fig. 1 Gamboa manufacturing software [Quintanilla, 2015] current implementation consists of two applications which exchange informations through sockets. Both applications can be associated either to an emulation on Arena or to the real workshop via its logic controller [Gamboa Quintanilla et al., 2016]. Both are intimately coupled with the actual workshop contents so that any workshop modification implies a software review and maintenance. The programs includes 160 classes, 1240 methods and 14802 lines of code.

A study of the software architecture seen of Figure 1 pointed out several limits from a software quality point of view, including major trends such as maintenance and evolution, verification, adaptability and reconfiguration, etc.

• Abstraction: Except some limited UML documentation, the application is mainly java code which means that (i) the evolutions must be performed by experimented software developers, (ii) each workshop modification implies a serious development and verification work, (iii) during verification, business problem issues are melt with implementation issues (iv) the general mechanisms are not distinct from those specific to the workshop and lead to no simple reuse, no capitalization, no simple reconfiguration... • Maintenance and Evolution: Agility is widespread now, users requires new requirements quickly. As the current manufacturing systems are highly diverse according to a mechanical or physical requirements, e.g., adding, modifying or deleting a mechanical component from the production workshop, any change can have a negative influence on the control software, which makes its maintenance and evolution a very delicate task. • Verification: As the complexity of control software for manufacturing systems shall increase in the next decade, proving the correctness of such properties before deploying the application in a real production cell becomes a critical task. The current framework does not offer a model in which we can apply the verification techniques and ensure the requirements fulfillment such as : correctness, QoS, consistency, etc. • Adaptability and Reuse Return on invest (ROI) applies not only when building workshops but also for the control software. The models and applications should be reused in different contexts. The current software architecture has not been implemented with an explicit and modular way that facilitates to reuse (parts of) the existing platform in different (rather similar) contexts.

The new vision intends to contribute to reduce these limitations.

Proposed approach

The proposal targets two goals: improve the quality of the software architecture (modularity, extendibility, variability) and improve the software construction process (abstraction, reuse, testability, reconfigurability).

Application Architecture

Instead of monolithic, nested and specific applications (left part of Figure 2) we target modular, independent and parametric applications (right part of Figure 2).

The goal is to be as modular as possible in order to reduce the coupling between the different components, and reduce the specificity criteria of a production software to a precise context (layout, orders, resources, etc.). The main and first issue is to separate the generic part, that stands for any kind of workshop, from the specific part which depends on the workshop resources and organization to put in place. The second issue is to separate the workshop features (resources and flows) from the manufacturing features (orders and products). The primer are rather static (even when considering reconfigurable manufacturing systems) while the latter are merely dynamic, time depending and subject to quality of service constraints. The GUI layer is separated in three applications: the system definition, configuration and monitoring. It is set up for constructing all interactions that an expert user can perform on a given manufacturing system, such as defining the various scenarios (Manufacturing editor), monitoring the running scenarios and define the layout on which each scenario will unfold. Technically, a web GUI layer enables not only to be independent from the operating systems and physical devices, but also a compliant property for connecting the applications to service-oriented cloud or IoT and Scada platforms. The running SoHMS application drives the production system under the given constraints of the manufacturing process. The process engine, that includes the HMS definition, scheduling algorithms and properties, is generic while the orders, products and key indicators are specific. Finally the physical layer corresponds to different needs: emulation with Arena for example, management blackboards and of course the real workshop. The communication network between the layers should be compliant with different implementations (sockets, messages, shared memories...) according to the target infrastructure.

Software Construction Process

Basically, the existing SoHMS application (left part of Figure 3) consists of a collection of Java programs and workshop description which have been manually written and tested for a specific workshop case study. The new vision (right part of Figure 3) installs a kind of software product line where people reason about the workshop to put in place before implementing it. HMS libraries are defined through components and services with related implementations, such that building manufacturing systems is rather an assembly and integra-tion activity rather than a programming activity. Handling models enables to work with abstract concepts and to verify the expected properties without taking into account useless implementation details. These are qualified as platform independent models (PIM) in the MDA spirit while the current programs are really platform specific models (PSM). Figure 3 illustrates the proposed approach using three layers:

• Modelling layer: Working with abstract models enables to address the problem complexity and to formalize the system's structure, behaviour and requirements. Following the principle of "separation of concerns", building orthogonal models enforces consistency and reduces the software complexity. As an example, we have separated the development of workshop models that specify the structure of the production cells and the paths supporting the product transportation from the manufacturing models that specify the product recipes, the orders management and the orchestration of the process. • Engineering layer: Model driven engineering means considering the model as a first class element in the development process. In this context, we aim to use the abstract models that we have defined in the modelling layer to apply : (i) formal verification, to get a high level of insurance about the correctness in early phase of the software development life cycle, (ii) models refinement, to automatically reduce the abstraction level until generate a concrete model that represents the application code, which increases the automation level in our approach. In addition we aim to construct standard and generic libraries and frameworks that can be integrated during the refinement process. • Implementation layer: Various implementation targets but also various goals can be motivated. The implementation layer will support the automation process to generate the logic code for each concern designed in the modelling stage (GUI code, Emulation testbed, Manager blackboard and reporting, etc.).

Abstraction and separation of concerns are key issues for quality software design.

Experimentation

This section will briefly introduce the application context, an illustrative example and will include various references on previous works and experimentation. A very basic software production chain version is implemented to illustrate the process of Figure 3. In this section we will present the resulted manufacturing software and the different choices of technologies to implement it.

Case Study and Methodology

We could restart from the case study of [Quintanilla, 2015] and refactor the existing programs. However the reverse engineering and refactoring process of this proof of concept (POC) application would take a long time because the specific and generic aspects of the solution are interleaved. This POC is still a resource for reverse-engineering but not the target system. So, the starting point for our case study was chosen as the Bench4Star [START_REF] Trentesaux | Benchmarking flexible job-shop scheduling and control systems[END_REF] manufacturing system. It is composed of seven workstations placed around a flexible transportation system. Workstations can have specific or concurrent capabilities. The production program consists of assembling components on a plate. From a methodological point of view, we follow an agile perspective and will focus on the transportation issues for the first iterations. In parallel to this engineering stage, a reverse engineering process is applied to the current [Quintanilla, 2015] application in order to extract the main components and separate the generic and specific abstractions. This work provides fruitful components to build the generic SoHMS framework of Figure 3.

Workshop and Manufacturing Editor

The manufacturing editor, located in the GUI layer of Figure 3, aims to describe the workshop features (resources, layouts, products and orders). In its first version, a simple web-based textual editor has been developed to define the manufacturing scenarios to play later on (Figure 4). Scenarios are defined using JSON data format, a widespread data exchange language, easy to read or write for humans and simple to be parsed by different programming languages. The text editor is implemented using different web technologies such as Vue.js, Bootstrap 4, HTML5, CodeMirror, WebScoket and PHP7. Details on this editor can be found in [START_REF] Carat | Contribution au logiciel de contrôle d'une flotte d'AGVs[END_REF]. The future versions will include GUI editors with menus and input forms equipped with drag and drop visual facilities to build workshops and manufacturing models.

Model Engineering (UML models)

For the modeling layer, the workshop and manufacturing models will be described using standard or domain specific languages (DSL) like UML, OCL, SysML, AADL, Kmelia as mentioned in [START_REF] André | Trusted Services for Cyber Manufacturing Systems[END_REF]. In the MDE spirit, the idea is to use a collection of languages, each of them enabling to capture a part of the system's semantics, say a point of view, that could serve as a support to reason about the system and to refine into implementable solutions.

Currently, we used UML diagrams to model elements of the manufacturing workshop and its behaviour. To find these elements, we got inspired by the application of [Quintanilla, 2015] and Bench4Star [START_REF] Trentesaux | Benchmarking flexible job-shop scheduling and control systems[END_REF]. For example, the class diagram of Figure 5 illustrates a simplified structural representation for the workshop (red part) and the production orders and scenarios (blue part). To specify the behaviour of the system, we used UML sequence diagrams to express the interactions and statecharts to represent the component dynamic evolution. Details on the current modeling activity can be found in [START_REF] Carat | Contribution au logiciel de contrôle d'une flotte d'AGVs[END_REF].

Model Transformation and Code Generation

The goal is to reduce the complexity of the development (programming and testing) by an iterative process made of simple processing task (modelling, verification, generation) i.e. small steps instead of big steps. Model transformation and code generation help to automate the refinement from high level abstract models to low level models (programs) that are specific to the target implementations. As an example, the JSON files of section 4.2 are used to generate the SoHMS specific elements. To generate the SoHMS code from the defined UML diagrams, we found that Papyrus1 , a graphical Modeling tool for UML with model transformation facilities, is a suitable tool to generate the corresponding code of our defined models basing on Model driven Engineering approach and using the model transformation concepts.

Toward Generic Frameworks and Libraries

We started a reverse-engineering activity that separate in [Quintanilla, 2015] application those classes that are generic, i.e. that are independent from the workshop to put in place, from those that are specific. The generic parts will be included in the SoHMS framework and the WebGUI framework: component libraries, HMS engine, default ordering policies, communication support (WebSocket technology), etc. The goal is to provide a generic and flexible communication support that can be extended to various platforms for developing such control systems.

Discussion

The use of MDE was already demonstrated on classical manufacturing control systems to improve productivity of over 50% compared to a hand-made code generation [START_REF] Cuadrado | Applying model-driven engineering in small software enterprises[END_REF], in flexibility and reconfigurability for SMEs [START_REF] Masood | A model-driven approach to enabling change capability in SMEs[END_REF] and in the application of products customization principles [START_REF] Aleksic | A case study on the object-oriented framework for modeling product families with the dominant variation of the topology in the one-of-a-kind production[END_REF]. On the other hand, numerous issues in innovative CPPS control systems can be pointed out, among which:

• Lack of performance analysis and guarantees, due to the emergent behaviour targeted by the application of multi-agent concepts. It is then necessary to implement trustable validation and verifications techniques [START_REF] Nastov | Mbse and v&v: a tool-equipped method for combining various v&v strategies[END_REF] in order to guarantee the respect of the design constraints in the final code; • Lack of a repeatable methodology, in order to stabilize the code generation process and increase the trust in the final behaviour of the system. This question is of major matter considering ethical questions, such as responsibility in case of major damages caused by the system [START_REF] Trentesaux | Designing ethical cyberphysical industrial systems[END_REF];

• Lack of stability of the underlying algorithms: in a Industry 4.0 context where service-oriented control systems are among the main objectives, the current developments need a major and constant rewriting of the code which prevent any reuse of codes or algorithms; • Difficulty in migrating from and managing legacy environments: the control system's evolution is of a great concern when implementing innovative control systems to existing workshops. Indeed, the investment costs makes it quite impossible to change the resources accordingly with the software. Therefore, the software needs to adapt without rewriting code to a frequent reconfiguration of resources along the evolutions of the workshop layout and the customers' demand.

Considering this context, the objectives of this proposal are in two ways:

• Building a prototype software production chain in the context of next generation manufacturing control systems using MDE as a POC. By doing so, generic models and libraries are to be developed and constructed that will constitute a step beyond current literature in terms of genericity of CPPS control systems development.

• Develop the advantages of the application of MDE in this POC, in order to be able to demonstrate the Validation and Verification possibilities of such an approach and its benefits on the global process of control systems development: this objective is directly bound to the richness of the models used for generation.

Conclusion and Perspectives

Introducing the service orientation into CPPS control software was a first step to build scalable, modular and integrated control software. The next step is to introduce model driven engineering to build modular software systems that are easier to reconfigure, to maintain, but also to reuse parts into new other systems. We introduce the new construction systems and experimented parts to show its feasibility. Future work will instrument the full software product line that implements our vision. At first, we will build the application engines i.e. the generic frameworks. Second, thanks to model transformation, we can generate the plugable software elements that are specific to the workshop to control. Third, we will write the code generators that enable to target different physical systems. Last, but not least, we will investigate the verification issues in order to implement trusty software systems.

Fig. 2

 2 Fig. 2 SoHMS Software Architecture

Fig. 3

 3 Fig. 3 Software Construction Architecture

Fig. 4

 4 Fig. 4 Web GUI first implementation for CPPS control software

Fig. 5

 5 Fig. 5 SoHMS UML Models

https://www.eclipse.org/papyrus/